text:
C***********************************************************************
C***********************************************************************
C
C   Version:        0.3
C   Last modified:  December 27, 1994
C   Authors:        Esmond G. Ng and Barry W. Peyton
C
C   Mathematical Sciences Section, Oak Ridge National Laboratory
C
C***********************************************************************
C***********************************************************************
C****************    FSUP1 ..... FIND SUPERNODES #1    *****************
C***********************************************************************
C***********************************************************************
C
C   PURPOSE:
C       THIS SUBROUTINE IS THE FIRST OF TWO ROUTINES FOR FINDING A
C       MAXIMAL SUPERNODE PARTITION.  IT RETURNS ONLY THE NUMBER OF
C       SUPERNODES NSUPER AND THE SUPERNODE MEMBERSHIP VECTOR SNODE(*),
C       WHICH IS OF LENGTH NEQNS.  THE VECTORS OF LENGTH NSUPER ARE
C       COMPUTED SUBSEQUENTLY BY THE COMPANION ROUTINE FSUP2.
C
C   METHOD AND ASSUMPTIONS:
C       THIS ROUTINE USES THE ELIMINATION TREE AND THE FACTOR COLUMN
C       COUNTS TO COMPUTE THE SUPERNODE PARTITION; IT ALSO ASSUMES A
C       POSTORDERING OF THE ELIMINATION TREE.
C
C   INPUT PARAMETERS:
C       (I) NEQNS       -   NUMBER OF EQUATIONS.
C       (I) ETPAR(*)    -   ARRAY OF LENGTH NEQNS, CONTAINING THE
C                           ELIMINATION TREE OF THE POSTORDERED MATRIX.
C       (I) COLCNT(*)   -   ARRAY OF LENGTH NEQNS, CONTAINING THE
C                           FACTOR COLUMN COUNTS: I.E., THE NUMBER OF
C                           NONZERO ENTRIES IN EACH COLUMN OF L
C                           (INCLUDING THE DIAGONAL ENTRY).
C
C   OUTPUT PARAMETERS:
C       (I) NOFSUB      -   NUMBER OF SUBSCRIPTS.
C       (I) NSUPER      -   NUMBER OF SUPERNODES (<= NEQNS).
C       (I) SNODE(*)    -   ARRAY OF LENGTH NEQNS FOR RECORDING
C                           SUPERNODE MEMBERSHIP.
C
C   FIRST CREATED ON    JANUARY 18, 1992.
C   LAST UPDATED ON     NOVEMBER 11, 1994.
C
C***********************************************************************
C
      SUBROUTINE  FSUP1  (  NEQNS , ETPAR , COLCNT, NOFSUB, NSUPER,
     &                      SNODE                                   )
C
C***********************************************************************
C
C       -----------
C       PARAMETERS.
C       -----------
        INTEGER             NEQNS , NOFSUB, NSUPER
        INTEGER             COLCNT(*)     , ETPAR(*)      ,
     &                      SNODE(*)
C
C       ----------------
C       LOCAL VARIABLES.
C       ----------------
        INTEGER             KCOL
C
C***********************************************************************
C
C       --------------------------------------------
C       COMPUTE THE FUNDAMENTAL SUPERNODE PARTITION.
C       --------------------------------------------
        NSUPER = 1
        SNODE(1) = 1
        NOFSUB = COLCNT(1)
        DO  300  KCOL = 2, NEQNS
            IF  ( ETPAR(KCOL-1) .EQ. KCOL )  THEN
                IF  ( COLCNT(KCOL-1) .EQ. COLCNT(KCOL)+1 )  THEN
                    SNODE(KCOL) = NSUPER
                    GO TO 300
                ENDIF
            ENDIF
            NSUPER = NSUPER + 1
            SNODE(KCOL) = NSUPER
            NOFSUB = NOFSUB + COLCNT(KCOL)
  300   CONTINUE
C
      RETURN
      END

max_stars_repo_name: Pangqiyuangh/SeIInv

max_issues_repo_name: Pangqiyuangh/SeIInv

max_forks_repo_name: Pangqiyuangh/SeIInv


text:
      FUNCTION RAN2(IDUM)
      PARAMETER (M=714025,IA=1366,IC=150889,RM=1.4005112E-6)
      DIMENSION IR(97)
      DATA IFF /0/
      IF(IDUM.LT.0.OR.IFF.EQ.0)THEN
        IFF=1
        IDUM=MOD(IC-IDUM,M)
        DO 11 J=1,97
          IDUM=MOD(IA*IDUM+IC,M)
          IR(J)=IDUM
11      CONTINUE
        IDUM=MOD(IA*IDUM+IC,M)
        IY=IDUM
      ENDIF
      J=1+(97*IY)/M
      IF(J.GT.97.OR.J.LT.1)PAUSE
      IY=IR(J)
      RAN2=IY*RM
      IDUM=MOD(IA*IDUM+IC,M)
      IR(J)=IDUM
      RETURN
      END

ext: for


text:
* $Id$
      SUBROUTINE txs_SECOND(T)
      common /intgop/ ncache,maxprice,iprint,iblock
      double precision T
c     REAL*8 T
      double precision util_cpusec
      external util_cpusec
*
*     rjh ... no point timing if not printing
*
      t = 0.d0
      if (iprint .le. 0) return
*
C      DIMENSION ITIME(4)
C     I1=TIMES(ITIME)
C     WRITE(*,*) ITIME
C     T=DBLE(ITIME(1)+ITIME(2)+ITIME(3)+ITIME(4))/100.0D0
C      T=MCLOCK()/100.0D0
      T = util_cpusec()
C
      RETURN
      END
C
*rak:      SUBROUTINE FDATE(X)
*rak:      CHARACTER*24 BLANK,X
*rak:      DATA BLANK /'                        '/
*rak:      X=BLANK
*rak:      END

ext: f

max_stars_repo_name: dinisAbranches/nwchem

max_issues_repo_name: dinisAbranches/nwchem

max_forks_repo_name: dinisAbranches/nwchem


text:
	subroutine indexx(n, arrin, indx)

	implicit none

c	Parameters:
	integer	n
	real	arrin(n)
	integer	indx(n)

c	Local variables:
	integer	i
	integer	indxt
	integer	ir
	integer	j
	integer	l
	real	q

      if (n.lt.1) then
         return
      else if (n.eq.1) then
         indx(1) = 1
         return
      endif

      DO 11 J=1,N
        INDX(J)=J
11    CONTINUE
      L=N/2+1
      IR=N
10    CONTINUE
        IF(L.GT.1)THEN
          L=L-1
          INDXT=INDX(L)
          Q=ARRIN(INDXT)
        ELSE
          INDXT=INDX(IR)
          Q=ARRIN(INDXT)
          INDX(IR)=INDX(1)
          IR=IR-1
          IF(IR.EQ.1)THEN
            INDX(1)=INDXT
            RETURN
          ENDIF
        ENDIF
        I=L
        J=L+L
20      IF(J.LE.IR)THEN
          IF(J.LT.IR)THEN
            IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1
          ENDIF
          IF(Q.LT.ARRIN(INDX(J)))THEN
            INDX(I)=INDX(J)
            I=J
            J=J+J
          ELSE
            J=IR+1
          ENDIF
        GO TO 20
        ENDIF
        INDX(I)=INDXT
      GO TO 10
      END


text:
      subroutine evdat4(arr,ipc,ipcxmx,k,nmax,narxmx,narxt,ntpr,
     $ ntprmx,prop,tempc)
c
c     This subroutine evaluates a thermodynamic property as a function
c     of temperature, using an interpolating polynomial whose
c     coefficients are stored in a 4D array arr. The second dimension
c     of this array corresponds to a temperature range. The third
c     dimension of this array usually corresponds to an order parameter,
c     as for pressure correction.
c
c     Compare with EQLIB/evdat3.f, in which arr is a 3D array, and
c     EQLIB/evdat3.f, in which this is a 3D array.
c
c     This subroutine is called by:
c
c       EQLIB/evdatr.f
c
c-----------------------------------------------------------------------
c
c     Principal input:
c
c       k      = the index for the fourth dimension of the arr array
c       arr    = two dimensional array of polynomial coefficients
c                  describing some thermodynamic function
c       ipcxmx = the third dimension of the arr array
c       ipc    = the index for the third dimension of the arr array
c       tempc  = temperature, C
c       ntpr   = temperature range flag
c       narxmx = first dimension of the arr array, the number
c                  of coefficients per temperature range
c       narxt  = array of numbers of coefficients in temperature
c                  ranges
c       ntprmx = second dimension of the arr array, the number
c                  of temperature ranges.
c       nmax   = third dimension of the arr array (can be ngtmax,
c                  nmtmax, or ngtmax)
c
c     Prinicpal output:
c
c       prop   = the calculated property array
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer ipcxmx,nmax,narxmx,ntprmx
c
      integer narxt(ntprmx)
c
      integer ipc,k,ntpr
c
      real*8 arr(narxmx,ntprmx,ipcxmx,nmax),prop
      real*8 tempc
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer n,nn,nt
c
c-----------------------------------------------------------------------
c
c     Evaluate the polynomial.
c
      prop = 0.
      nt = narxt(ntpr)
      do nn = 1,nt
        n = nt + 1 - nn
        prop = arr(n,ntpr,ipc,k) + tempc*prop
      enddo
c
      end


text:
C--- SPARSPAK-A (ANSI FORTRAN) RELEASE III --- NAME = FNROOT
C  (C)  UNIVERSITY OF WATERLOO   JANUARY 1984
C***************************************************************
C***************************************************************
C******     FNROOT ..... FIND PSEUDO-PERIPHERAL NODE     *******
C***************************************************************
C***************************************************************
C
C     PURPOSE - FNROOT IMPLEMENTS A MODIFIED VERSION OF THE
C        SCHEME BY GIBBS, POOLE, AND STOCKMEYER TO FIND PSEUDO-
C        PERIPHERAL NODES.  IT DETERMINES SUCH A NODE FOR THE
C        SECTION SUBGRAPH SPECIFIED BY MASK AND ROOT.
C
C     INPUT PARAMETERS -
C        (XADJ,ADJNCY) - ADJACENCY STRUCTURE PAIR FOR THE GRAPH.
C        MASK   - SPECIFIES A SECTION SUBGRAPH. NODES FOR WHICH
C                 MASK IS ZERO ARE IGNORED BY FNROOT.
C
C     UPDATED PARAMETER -
C        ROOT   - ON INPUT, IT (ALONG WITH MASK) DEFINES THE
C                 COMPONENT FOR WHICH A PSEUDO-PERIPHERAL NODE
C                 IS TO BE FOUND.  ON OUTPUT, IT IS THE NODE
C                 OBTAINED.
C
C     OUTPUT PARAMETERS -
C        NLVL   - IS THE NUMBER OF LEVELS IN THE LEVEL
C                 STRUCTURE ROOTED AT THE NODE ROOT.
C        (XLS,LS) - THE LEVEL STRUCTURE ARRAY PAIR CONTAINING
C                 THE LEVEL STRUCTURE FOUND.
C
C    PROGRAM SUBROUTINE -
C       ROOTLS.
C
C***************************************************************
C
      SUBROUTINE  FNROOT ( ROOT, XADJ, ADJNCY, MASK, NLVL,
     1                     XLS, LS )
C
C***************************************************************
C
         INTEGER    ADJNCY(1), LS(1)    , MASK(1)  , XLS(1)
         INTEGER    XADJ(1)
         INTEGER    CCSIZE, J     , JSTRT , K     , KSTOP ,
     1              KSTRT , MINDEG, NABOR , NDEG  , NLVL  ,
     1              NODE  , NUNLVL, ROOT
C
C***************************************************************
C
C        ---------------------------------------------
C        DETERMINE THE LEVEL STRUCTURE ROOTED AT ROOT.
C        ---------------------------------------------
         CALL  ROOTLS ( ROOT, XADJ, ADJNCY, MASK, NLVL, XLS, LS )
         CCSIZE = XLS(NLVL+1) - 1
         IF  ( NLVL .EQ. 1  .OR.  NLVL .EQ. CCSIZE )  RETURN
C        ----------------------------------------------------
C        PICK A NODE WITH MINIMUM DEGREE FROM THE LAST LEVEL.
C        ----------------------------------------------------
  100    CONTINUE
             JSTRT = XLS(NLVL)
             MINDEG = CCSIZE
             ROOT = LS(JSTRT)
             IF  ( CCSIZE .EQ. JSTRT )  GO TO 400
                 DO  300  J = JSTRT, CCSIZE
                     NODE = LS(J)
                     NDEG = 0
                     KSTRT = XADJ(NODE)
                     KSTOP = XADJ(NODE+1) - 1
                     DO  200  K = KSTRT, KSTOP
                         NABOR = ADJNCY(K)
                         IF  ( MASK(NABOR) .GT. 0 )  NDEG = NDEG + 1
  200                CONTINUE
                     IF  ( NDEG .GE. MINDEG )  GO TO 300
                         ROOT = NODE
                         MINDEG = NDEG
  300            CONTINUE
  400        CONTINUE
C            ----------------------------------------
C            AND GENERATE ITS ROOTED LEVEL STRUCTURE.
C            ----------------------------------------
             CALL  ROOTLS ( ROOT, XADJ, ADJNCY, MASK,
     1                      NUNLVL, XLS, LS )
             IF  ( NUNLVL .LE. NLVL )  RETURN
             NLVL = NUNLVL
             IF  ( NLVL .LT. CCSIZE )  GO TO 100
                 RETURN
C
      END


program main
    use FEMSolverClass
    implicit none

    type(FEMSolver_) :: solver
    type(FEMDomain_),target :: domains(1)
    type(IO_) :: f
    integer(int32) :: i
    real(real64) :: Vs,t,dt,E_Al
    real(real64),allocatable :: Mode_U(:),mode_Ut(:),freq(:),eigen_value(:),eigen_vectors(:,:)

    ! Modal analysis
!Create file or...
    call domains(1)%create("Cube3D",x_num=2,y_num=6,z_num=40)
    call domains(1)%resize(x=1.0d0)
    call domains(1)%resize(y=5.0d0)
    call domains(1)%resize(z=50.0d0)
    call domains(1)%vtk("AlminiumBar")

!read file
    !call domains(1) % read("AlminiumBar.vtk")

    call solver%init(NumDomain=1,NumInterfaceElement=0)
    call solver%setDomain(FEMDomain=domains(1),DomainID=1)
    call solver%setCRS(DOF=3)

    !$OMP parallel do
    do i=1,domains(1)%ne()
        call solver%setMatrix(&
            DomainID=1,&
            ElementID=i,&
            DOF=3,&
            Matrix=domains(1)%MassMatrix(ElementID=i,Density=2.70d0,DOF=3) &
            )
    enddo
    !$OMP end parallel do

    call solver%keepThisMatrixAs("B")
    call solver%zeros()

    print *, "Save Stiffness Matrix"

    E_Al = 70.0d0*1000.0d0*1000.0d0
    !$OMP parallel do
    do i=1,domains(1)%ne()
        call solver%setMatrix(&
            DomainID=1,&
            ElementID=i,&
            DOF=3,&
            Matrix=domains(1)%StiffnessMatrix(ElementID=i,E=70.0d0*1000.0d0*1000.0d0,v=0.300d0) &
            )
    enddo
    !$OMP end parallel do

    call solver%keepThisMatrixAs("A")

    ! Eigen value problem solver by scipy

    print *, "solver%eig"
    call solver%eig(num_eigen=domains(1)%nn()*domains(1)%nd(),eigen_value=eigen_value,eigen_vectors=eigen_vectors)

    ! read results
    freq = sqrt(abs(eigen_value))/2.0d0/3.141590d0
    dt = 0.10d0

    ! 20 modes
    do i_i=1,20
        mode_U = zeros(size(eigen_vectors,1))
        mode_U = eigen_vectors(:,i_i)
        dt = 1.0d0/freq(i_i)/100.0d0
        do j_j=1,100
            t = dt * dble(j_j-1)
            mode_Ut = mode_U*cos( 2.0d0*3.140d0*freq(i_i)*t )

            domains(1)%mesh%nodcoord = domains(1)%mesh%nodcoord &
            +10.0d0*reshape(mode_Ut,domains(1)%nn(),3 )

            call domains(1)%vtk("Mode_Fortran_"+str(i_i)+"_t_"+str(j_j))
            domains(1)%mesh%nodcoord = domains(1)%mesh%nodcoord &
            -10.0d0*reshape(mode_Ut,domains(1)%nn(),3 )
        enddo
    enddo

    print *, "Freq (Hz)"
    call print(sqrt(abs(eigen_value(1:30)))/2.0d0/3.141590d0 )

end program main


text:
module Mod_AdaptiveHexa
   use Mod_AdaptiveInterpolationList
   implicit none

contains

   subroutine TrilinearHexa_GetElementVariation(elcod,ivariation)
      use typre
      implicit none
      real(rp) :: elcod(:,:)
      integer(ip) :: ivariation

      ivariation = 1

   end subroutine

   subroutine TrilinearHexa_ElementGetDimensions(pnodb,nchild,nface,nchildface)
      use typre
      implicit none
      integer(ip), intent(out) :: pnodb,nchild,nface,nchildface

      pnodb = 4
      nchild = 8
      nface = 6
      nchildface = 4
   end subroutine

   subroutine TrilinearHexa_Add_node_to_child(nchild,inode,jnode,pnods,lnods,ipoin,PointerToChildren)
      !This subroutine adds a new node to a child element when this is built
      use typre
      implicit none
      integer(ip), intent(in) :: nchild,inode,jnode,ipoin
      integer(ip), intent(inout) :: pnods(*),lnods(*),PointerToChildren(*)

      integer(ip) :: knode
      integer(ip) :: ielem,ispos

      if (inode == jnode) then
         knode = inode
      elseif ((inode == 1 .and. jnode == 2) .or. (inode == 2 .and. jnode == 1)) then
         knode = 9
      elseif ((inode == 1 .and. jnode == 4) .or. (inode == 4 .and. jnode == 1)) then
         knode = 10
      elseif ((inode == 1 .and. jnode == 5) .or. (inode == 5 .and. jnode == 1)) then
         knode = 11
      elseif ((inode == 2 .and. jnode == 3) .or. (inode == 3 .and. jnode == 2)) then
         knode = 12
      elseif ((inode == 2 .and. jnode == 6) .or. (inode == 6 .and. jnode == 2)) then
         knode = 13
      elseif ((inode == 3 .and. jnode == 4) .or. (inode == 4 .and. jnode == 3)) then
         knode = 14
      elseif ((inode == 3 .and. jnode == 7) .or. (inode == 7 .and. jnode == 3)) then
         knode = 15
      elseif ((inode == 4 .and. jnode == 8) .or. (inode == 8 .and. jnode == 4)) then
         knode = 16
      elseif ((inode == 5 .and. jnode == 6) .or. (inode == 6 .and. jnode == 5)) then
         knode = 17
      elseif ((inode == 5 .and. jnode == 8) .or. (inode == 8 .and. jnode == 5)) then
         knode = 18
      elseif ((inode == 6 .and. jnode == 7) .or. (inode == 7 .and. jnode == 6)) then
         knode = 19
      elseif ((inode == 7 .and. jnode == 8) .or. (inode == 8 .and. jnode == 7)) then
         knode = 20
      endif

      select case (knode)
      case(1)
         ielem = PointerToChildren(1)
         ispos = pnods(ielem)-1
         lnods(ispos+1) = ipoin
      case(2)
         ielem = PointerToChildren(5)
         ispos = pnods(ielem)-1
         lnods(ispos+2) = ipoin
      case(3)
         ielem = PointerToChildren(7)
         ispos = pnods(ielem)-1
         lnods(ispos+3) = ipoin
      case(4)
         ielem = PointerToChildren(3)
         ispos = pnods(ielem)-1
         lnods(ispos+4) = ipoin
      case(5)
         ielem = PointerToChildren(2)
         ispos = pnods(ielem)-1
         lnods(ispos+5) = ipoin
      case(6)
         ielem = PointerToChildren(6)
         ispos = pnods(ielem)-1
         lnods(ispos+6) = ipoin
      case(7)
         ielem = PointerToChildren(8)
         ispos = pnods(ielem)-1
         lnods(ispos+7) = ipoin
      case(8)
         ielem = PointerToChildren(4)
         ispos = pnods(ielem)-1
         lnods(ispos+8) = ipoin
      case(9)
         ielem = PointerToChildren(1)
         ispos = pnods(ielem)-1
         lnods(ispos+2) = ipoin

         ielem = PointerToChildren(5)
         ispos = pnods(ielem)-1
         lnods(ispos+1) = ipoin
      case(10)
         ielem = PointerToChildren(1)
         ispos = pnods(ielem)-1
         lnods(ispos+4) = ipoin

         ielem = PointerToChildren(3)
         ispos = pnods(ielem)-1
         lnods(ispos+1) = ipoin
      case(11)
         ielem = PointerToChildren(1)
         ispos = pnods(ielem)-1
         lnods(ispos+5) = ipoin

         ielem = PointerToChildren(2)
         ispos = pnods(ielem)-1
         lnods(ispos+1) = ipoin
      case(12)
         ielem = PointerToChildren(5)
         ispos = pnods(ielem)-1
         lnods(ispos+3) = ipoin

         ielem = PointerToChildren(7)
         ispos = pnods(ielem)-1
         lnods(ispos+2) = ipoin
      case(13)
         ielem = PointerToChildren(5)
         ispos = pnods(ielem)-1
         lnods(ispos+6) = ipoin

         ielem = PointerToChildren(6)
         ispos = pnods(ielem)-1
         lnods(ispos+2) = ipoin
      case(14)
         ielem = PointerToChildren(3)
         ispos = pnods(ielem)-1
         lnods(ispos+3) = ipoin

         ielem = PointerToChildren(7)
         ispos = pnods(ielem)-1
         lnods(ispos+4) = ipoin
      case(15)
         ielem = PointerToChildren(7)
         ispos = pnods(ielem)-1
         lnods(ispos+7) = ipoin

         ielem = PointerToChildren(8)
         ispos = pnods(ielem)-1
         lnods(ispos+3) = ipoin
      case(16)
         ielem = PointerToChildren(3)
         ispos = pnods(ielem)-1
         lnods(ispos+8) = ipoin

         ielem = PointerToChildren(4)
         ispos = pnods(ielem)-1
         lnods(ispos+4) = ipoin
      case(17)
         ielem = PointerToChildren(2)
         ispos = pnods(ielem)-1
         lnods(ispos+6) = ipoin

         ielem = PointerToChildren(6)
         ispos = pnods(ielem)-1
         lnods(ispos+5) = ipoin
      case(18)
         ielem = PointerToChildren(2)
         ispos = pnods(ielem)-1
         lnods(ispos+8) = ipoin

         ielem = PointerToChildren(4)
         ispos = pnods(ielem)-1
         lnods(ispos+5) = ipoin
      case(19)
         ielem = PointerToChildren(6)
         ispos = pnods(ielem)-1
         lnods(ispos+7) = ipoin

         ielem = PointerToChildren(8)
         ispos = pnods(ielem)-1
         lnods(ispos+6) = ipoin
      case(20)
         ielem = PointerToChildren(4)
         ispos = pnods(ielem)-1
         lnods(ispos+7) = ipoin

         ielem = PointerToChildren(8)
         ispos = pnods(ielem)-1
         lnods(ispos+8) = ipoin
      end select

   end subroutine

   subroutine TrilinearHexa_get_edges_from_child(pnods,lnods,nchild,PointerToChildren,nedge,ledger)
      use typre
      implicit none
      integer(ip), intent(in) :: pnods(*),lnods(*),nchild,PointerToChildren(nchild)
      integer(ip), intent(out) :: nedge,ledger(3,*)

      integer(ip) :: ispos(8)

      nedge = 12

      ispos(:) = pnods(PointerToChildren(:))-1

      !first edge (node 1 2 9)
      ledger(1,1) = lnods(ispos(1)+1)
      ledger(2,1) = lnods(ispos(5)+2)
      ledger(3,1) = lnods(ispos(1)+2)

      !second edge (node 1 4 10)
      ledger(1,2) = lnods(ispos(1)+1)
      ledger(2,2) = lnods(ispos(3)+4)
      ledger(3,2) = lnods(ispos(1)+4)

      !third edge (1 5 11)
      ledger(1,3) = lnods(ispos(1)+1)
      ledger(2,3) = lnods(ispos(2)+5)
      ledger(3,3) = lnods(ispos(1)+5)

      !fourth edge (2 3 12)
      ledger(1,4) = lnods(ispos(5)+2)
      ledger(2,4) = lnods(ispos(7)+3)
      ledger(3,4) = lnods(ispos(5)+3)

      !fifth edge (2 6 13)
      ledger(1,5) = lnods(ispos(5)+2)
      ledger(2,5) = lnods(ispos(6)+6)
      ledger(3,5) = lnods(ispos(5)+6)

      !6 edge (3 4 14
      ledger(1,6) = lnods(ispos(7)+3)
      ledger(2,6) = lnods(ispos(3)+4)
      ledger(3,6) = lnods(ispos(7)+4)

      !7 edge (3 7 15)
      ledger(1,7) = lnods(ispos(7)+3)
      ledger(2,7) = lnods(ispos(8)+7)
      ledger(3,7) = lnods(ispos(7)+7)

      !8 (4 8 16)
      ledger(1,8) = lnods(ispos(3)+4)
      ledger(2,8) = lnods(ispos(4)+8)
      ledger(3,8) = lnods(ispos(3)+8)

      !9 (5 6 17)
      ledger(1,9) = lnods(ispos(2)+5)
      ledger(2,9) = lnods(ispos(6)+6)
      ledger(3,9) = lnods(ispos(2)+6)

      !10 (5 8 18)
      ledger(1,10) = lnods(ispos(2)+5)
      ledger(2,10) = lnods(ispos(4)+8)
      ledger(3,10) = lnods(ispos(2)+8)

      !11 edge (6 7 19)
      ledger(1,11) = lnods(ispos(6)+6)
      ledger(2,11) = lnods(ispos(8)+7)
      ledger(3,11) = lnods(ispos(6)+7)

      !12 edge (7 8 20)
      ledger(1,12) = lnods(ispos(8)+7)
      ledger(2,12) = lnods(ispos(4)+8)
      ledger(3,12) = lnods(ispos(8)+8)

   end subroutine

   subroutine TrilinearHexa_face_from_el(nface,pnodb,lnods,iface,sort,facnod)
      use typre
      implicit none

      integer(ip) :: iface,nface,pnodb
      integer(ip) :: facnod(pnodb,nface),lnods(*)
      integer(ip) :: sort

      integer(ip) :: faces(4,6), auxfaces(4),jface

      faces = reshape((/ 1,4,8,5, &
                         2,6,7,3, &
                         1,5,6,2, &
                         4,3,7,8, &
                         1,2,4,3, &
                         5,8,7,6  &
                     /),shape(faces))

      if (iface == 0) then
         do jface = 1,6
            auxfaces = faces(:,jface)
            facnod(:,jface) = lnods(auxfaces)
         enddo
         if (sort == 1) then
            do jface = 1,6
               call sortb(4,facnod(1,jface))
            enddo
         endif
      else
         auxfaces = faces(:,iface)
         facnod(:,1) = lnods(auxfaces)
         if (sort == 1) call sortb(4,facnod(1,1))
      endif

   end subroutine

   subroutine TrilinearHexa_up_faces(nface,ison,jface)
      !this subroutine gives the new neighbour face in the parent when an element is gone
      use typre
      implicit none
      integer(ip):: ielty,nface,jface,ison
      select case (ison)
         case (1)
            if (jface == 1) then
               jface = 1
            elseif (jface == 3) then
               jface = 3
            elseif (jface == 5) then
               jface = 5
            endif
         case (2)
            if (jface == 1) then
               jface = 1
            elseif (jface == 3) then
               jface = 3
            elseif (jface == 6) then
               jface = 6
            endif
         case (3)
            if (jface == 1) then
               jface = 1
            elseif (jface == 4) then
               jface = 4
            elseif (jface == 5) then
               jface = 5
            endif
         case (4)
            if (jface == 1) then
               jface = 1
            elseif (jface == 4) then
               jface = 4
            elseif (jface == 6) then
               jface = 6
            endif
         case (5)
            if (jface == 2) then
               jface = 2
            elseif (jface == 3) then
               jface = 3
            elseif (jface == 5) then
               jface = 5
            endif
         case (6)
            if (jface == 2) then
               jface = 2
            elseif (jface == 3) then
               jface = 3
            elseif (jface == 6) then
               jface = 6
            endif
         case (7)
            if (jface == 2) then
               jface = 2
            elseif (jface == 4) then
               jface = 4
            elseif (jface == 5) then
               jface = 5
            endif
         case (8)
            if (jface == 2) then
               jface = 2
            elseif (jface == 4) then
               jface = 4
            elseif (jface == 6) then
               jface = 6
            endif
         end select

   end subroutine

   subroutine TrilinearHexa_down_faces(jface,facchild)
   use typre
   implicit none
   integer(ip):: jface,facchild(2,*)
   !this subroutine gives the neighbour face in the child for a given face in the parent
      select case (jface)
      case (1)
         facchild(1,1) = 1    !first the son, then the face
         facchild(2,1) = 1
         facchild(1,2) = 2
         facchild(2,2) = 1
         facchild(1,3) = 3
         facchild(2,3) = 1
         facchild(1,4) = 4
         facchild(2,4) = 1
      case (2)
         facchild(1,1) = 5
         facchild(2,1) = 2
         facchild(1,2) = 6
         facchild(2,2) = 2
         facchild(1,3) = 7
         facchild(2,3) = 2
         facchild(1,4) = 8
         facchild(2,4) = 2
      case (3)
         facchild(1,1) = 1
         facchild(2,1) = 3
         facchild(1,2) = 2
         facchild(2,2) = 3
         facchild(1,3) = 5
         facchild(2,3) = 3
         facchild(1,4) = 6
         facchild(2,4) = 3
      case (4)
         facchild(1,1) = 3
         facchild(2,1) = 4
         facchild(1,2) = 4
         facchild(2,2) = 4
         facchild(1,3) = 7
         facchild(2,3) = 4
         facchild(1,4) = 8
         facchild(2,4) = 4
      case (5)
         facchild(1,1) = 1
         facchild(2,1) = 5
         facchild(1,2) = 3
         facchild(2,2) = 5
         facchild(1,3) = 5
         facchild(2,3) = 5
         facchild(1,4) = 7
         facchild(2,4) = 5
      case (6)
         facchild(1,1) = 2
         facchild(2,1) = 6
         facchild(1,2) = 4
         facchild(2,2) = 6
         facchild(1,3) = 6
         facchild(2,3) = 6
         facchild(1,4) = 8
         facchild(2,4) = 6
      end select

   end subroutine

   subroutine TrilinearHexa_nodes_to_face(pnodb,lnodb,iface)
      use typre
      implicit none
      integer(ip) :: ndime,pnode,pnodb,lnodb(pnodb),iface
      integer(ip) :: auxlnodb(4)

      auxlnodb = lnodb
      call sortb(4,auxlnodb)

      if (auxlnodb(1) == 1 .and. auxlnodb(2) == 4 .and. auxlnodb(3) == 5 .and. auxlnodb(4) == 8) then
         iface = 1
      elseif (auxlnodb(1) == 2 .and. auxlnodb(2) == 3 .and. auxlnodb(3) == 6 .and. auxlnodb(4) == 7) then
         iface = 2
      elseif (auxlnodb(1) == 1 .and. auxlnodb(2) == 2.and. auxlnodb(3) == 5 .and. auxlnodb(4) == 6) then
         iface = 3
      elseif (auxlnodb(1) == 3 .and. auxlnodb(2) == 4 .and. auxlnodb(3) == 7 .and. auxlnodb(4) == 8) then
         iface = 4
      elseif (auxlnodb(1) == 1 .and. auxlnodb(2) == 2 .and. auxlnodb(3) == 3 .and. auxlnodb(4) == 4) then
         iface = 5
      elseif (auxlnodb(1) == 5 .and. auxlnodb(2) == 6 .and. auxlnodb(3) == 7 .and. auxlnodb(4) == 8) then
         iface = 6
      endif

   end subroutine

   subroutine TrilinearHexa_ConnectInteriorFaces(jfpos,son,lfacs)
      use typre
      implicit none
      integer(ip) :: jfpos(*), son(*), lfacs(2,*)

      !interior faces
      !Elements 1 and 5
      lfacs(1,jfpos(1)+2) = son(5)
      lfacs(2,jfpos(1)+2) = 1

      lfacs(1,jfpos(5)+1) = son(1)
      lfacs(2,jfpos(5)+1) = 2

      !Elements 1 and 2
      lfacs(1,jfpos(1)+6) = son(2)
      lfacs(2,jfpos(1)+6) = 5

      lfacs(1,jfpos(2)+5) = son(1)
      lfacs(2,jfpos(2)+5) = 6

      !Elements  1 and 3
      lfacs(1,jfpos(1)+4) = son(3)
      lfacs(2,jfpos(1)+4) = 3

      lfacs(1,jfpos(3)+3) = son(1)
      lfacs(2,jfpos(3)+3) = 4

      !Elements 2 and 4
      lfacs(1,jfpos(2)+4) = son(4)
      lfacs(2,jfpos(2)+4) = 3

      lfacs(1,jfpos(4)+3) = son(2)
      lfacs(2,jfpos(4)+3) = 4

      !Elements 2 and 6
      lfacs(1,jfpos(2)+2) = son(6)
      lfacs(2,jfpos(2)+2) = 1

      lfacs(1,jfpos(6)+1) = son(2)
      lfacs(2,jfpos(6)+1) = 2

      !Elements 3 and 4
      lfacs(1,jfpos(3)+6) = son(4)
      lfacs(2,jfpos(3)+6) = 5

      lfacs(1,jfpos(4)+5) = son(3)
      lfacs(2,jfpos(4)+5) = 6

      !Elements 3 and 7
      lfacs(1,jfpos(3)+2) = son(7)
      lfacs(2,jfpos(3)+2) = 1

      lfacs(1,jfpos(7)+1) = son(3)
      lfacs(2,jfpos(7)+1) = 2

      !Elements 4 and 8
      lfacs(1,jfpos(4)+2) = son(8)
      lfacs(2,jfpos(4)+2) = 1

      lfacs(1,jfpos(8)+1) = son(4)
      lfacs(2,jfpos(8)+1) = 2

      !Elements 5 and 6
      lfacs(1,jfpos(5)+6) = son(6)
      lfacs(2,jfpos(5)+6) = 5

      lfacs(1,jfpos(6)+5) = son(5)
      lfacs(2,jfpos(6)+5) = 6

      !Elements 5 and 7
      lfacs(1,jfpos(5)+4) = son(7)
      lfacs(2,jfpos(5)+4) = 3

      lfacs(1,jfpos(7)+3) = son(5)
      lfacs(2,jfpos(7)+3) = 4

      !Elements 6 and 8
      lfacs(1,jfpos(6)+4) = son(8)
      lfacs(2,jfpos(6)+4) = 3

      lfacs(1,jfpos(8)+3) = son(6)
      lfacs(2,jfpos(8)+3) = 4

      !Elements 7 and 8
      lfacs(1,jfpos(7)+6) = son(8)
      lfacs(2,jfpos(7)+6) = 5

      lfacs(1,jfpos(8)+5) = son(7)
      lfacs(2,jfpos(8)+5) = 6

   end subroutine

   subroutine TrilinearHexa_AdaptiveGetChildrenExternalFacesList(inchild,jjchild,jjnface,jjface)
      use typre
      implicit none
      integer(ip) :: inchild,jjchild(:), jjnface(:), jjface(:,:)

      !there are 3 children with faces shared with the parent element
      inchild = 8
      !These elements are child 1, 2 and 3
      jjchild(1) = 1
      jjchild(2) = 2
      jjchild(3) = 3
      jjchild(4) = 4
      jjchild(5) = 5
      jjchild(6) = 6
      jjchild(7) = 7
      jjchild(8) = 8
      !Each child has 2 external faces
      jjnface(1:8)  = 3
      !The external faces are the 2 first ones for all the children
      jjface(1,1) = 1
      jjface(2,1) = 3
      jjface(3,1) = 5

      jjface(1,2) = 1
      jjface(2,2) = 3
      jjface(3,2) = 6

      jjface(1,3) = 1
      jjface(2,3) = 4
      jjface(3,3) = 5

      jjface(1,4) = 1
      jjface(2,4) = 4
      jjface(3,4) = 6

      jjface(1,5) = 2
      jjface(2,5) = 3
      jjface(3,5) = 5

      jjface(1,6) = 2
      jjface(2,6) = 3
      jjface(3,6) = 6

      jjface(1,7) = 2
      jjface(2,7) = 4
      jjface(3,7) = 5

      jjface(1,8) = 2
      jjface(2,8) = 4
      jjface(3,8) = 6
   end subroutine

   subroutine TrilinearHexa_GetInterpolationList(ivariation,InterpList)
      use typre
      implicit none
      integer(ip) :: ivariation
      type(InterpolationList) :: InterpList

      integer(ip) :: auxlnods(8)    , facnods(4),ichild,iface,inode

      auxlnods = (/1_ip, 2_ip, 3_ip, 4_ip, 5_ip, 6_ip, 7_ip, 8_ip /)

      InterpList%n = 7

      !The in the faces ones
      iface = 1
      ichild = 1
      inode  = 8
      call AddFace

      iface = 2
      ichild = 5
      inode  = 7
      call AddFace

      iface = 3
      ichild = 1
      inode  = 6
      call AddFace

      iface = 4
      ichild = 7
      inode  = 8
      call AddFace

      iface = 5
      ichild = 1
      inode  = 3
      call AddFace

      iface = 6
      ichild = 6
      inode  = 8
      call AddFace

      !The interior one
      InterpList%childElementList(7) = 1
      InterpList%childObjectiveNodeList(7) = 7
      InterpList%nInterpolatingNodes(7) = 8
      InterpList%InterpolatingNodesList(1:8,7) = (/1_ip, 2_ip, 3_ip, 4_ip, 5_ip, 6_ip, 7_ip, 8_ip /)
      InterpList%InterpolatingCoefficientsList(1:8,7) = 0.125_rp

contains
      subroutine AddFace
         implicit none

         InterpList%childElementList(iface) =   ichild
         InterpList%childObjectiveNodeList(iface) = inode
         InterpList%nInterpolatingNodes(iface) = 4
         call TrilinearHexa_face_from_el(6_ip,4_ip,auxlnods,iface,0,InterpList%InterpolatingNodesList(1:4,iface))
         InterpList%InterpolatingCoefficientsList(1:4,iface) = 0.25_rp
      end subroutine

   end subroutine

   subroutine TrilinearHexa_GetEdgeNodesForInode(inode,nnode,nodelist)
      use typre
      implicit none
      integer(ip) :: inode, nnode,nodelist(3)

      nnode = 3

      if (inode == 1) then
         nodelist(1) = 2
         nodelist(2) = 4
         nodelist(3) = 5
      elseif (inode == 2) then
         nodelist(1) = 1
         nodelist(2) = 3
         nodelist(3) = 6
      elseif (inode == 3) then
         nodelist(1) = 2
         nodelist(2) = 4
         nodelist(3) = 7
      elseif (inode == 4) then
         nodelist(1) = 1
         nodelist(2) = 3
         nodelist(3) = 8
      elseif (inode == 5) then
         nodelist(1) = 1
         nodelist(2) = 6
         nodelist(3) = 8
      elseif (inode == 6) then
         nodelist(1) = 2
         nodelist(2) = 5
         nodelist(3) = 7
      elseif (inode == 7) then
         nodelist(1) = 3
         nodelist(2) = 6
         nodelist(3) = 8
      elseif (inode == 8) then
         nodelist(1) = 4
         nodelist(2) = 5
         nodelist(3) = 7
      endif
   end subroutine

   subroutine TrilinearHexa_GetNumberOfInteriorNodesToAdd(nInteriorNodes)
      use typre
      implicit none
      integer(ip) :: nInteriorNodes

      nInteriorNodes = 1

   end subroutine

   subroutine TrilinearHexa_GetInteriorPointFromParent(inode,jchild,jnode)
      use typre
      implicit none
      integer(ip) :: inode,jchild,jnode

      if (inode == 1) then
         jchild = 1
         jnode = 7
      else
         call runend('Wrong number of interior nodes')
      endif
   end subroutine

   subroutine TrilinearHexa_AddNodeToChildInterior(ivariation,iadd,newpo,PointerToChildren,pnods,lnods)
      use typre
      implicit none
      integer(ip) :: ivariation,iadd,newpo,PointerToChildren(*),pnods(*),lnods(*)
      integer(ip) :: ispos(8), ichild,inode

      !Only one interior node needs to be added
      if (iadd == 1) then
         ichild = 1
         inode = 7
         ispos(ichild) = pnods(PointerToChildren(ichild))-1
         lnods(ispos(ichild)+inode) = newpo

         ichild = 2
         inode = 3
         ispos(ichild) = pnods(PointerToChildren(ichild))-1
         lnods(ispos(ichild)+inode) = newpo

         ichild = 3
         inode = 6
         ispos(ichild) = pnods(PointerToChildren(ichild))-1
         lnods(ispos(ichild)+inode) = newpo

         ichild = 4
         inode = 2
         ispos(ichild) = pnods(PointerToChildren(ichild))-1
         lnods(ispos(ichild)+inode) = newpo

         ichild = 5
         inode = 8
         ispos(ichild) = pnods(PointerToChildren(ichild))-1
         lnods(ispos(ichild)+inode) = newpo

         ichild = 6
         inode = 4
         ispos(ichild) = pnods(PointerToChildren(ichild))-1
         lnods(ispos(ichild)+inode) = newpo

         ichild = 7
         inode = 5
         ispos(ichild) = pnods(PointerToChildren(ichild))-1
         lnods(ispos(ichild)+inode) = newpo

         ichild = 8
         inode = 1
         ispos(ichild) = pnods(PointerToChildren(ichild))-1
         lnods(ispos(ichild)+inode) = newpo
      endif
   end subroutine

   subroutine TrilinearHexa_GetNedge(nedge)
      use typre
      implicit none
      integer(ip) :: nedge

      nedge = 12
   end subroutine

   subroutine TrilinearHexa_GetEdges(iedge,lnodedge)
      use typre
      implicit none
      integer(ip) :: ielty,iedge,lnodedge(*)

      select case(iedge)
         case(1)
            lnodedge(1) = 1
            lnodedge(2) = 2
         case(2)
            lnodedge(1) = 1
            lnodedge(2) = 4
         case(3)
            lnodedge(1) = 1
            lnodedge(2) = 5
         case(4)
            lnodedge(1) = 2
            lnodedge(2) = 3
         case(5)
            lnodedge(1) = 2
            lnodedge(2) = 6
         case(6)
            lnodedge(1) = 3
            lnodedge(2) = 4
         case(7)
            lnodedge(1) = 3
            lnodedge(2) = 7
         case(8)
            lnodedge(1) = 4
            lnodedge(2) = 8
         case(9)
            lnodedge(1) = 5
            lnodedge(2) = 6
         case(10)
            lnodedge(1) = 5
            lnodedge(2) = 8
         case(11)
            lnodedge(1) = 6
            lnodedge(2) = 7
         case(12)
            lnodedge(1) = 7
            lnodedge(2) = 8

      end select

   end subroutine

   subroutine TrilinearHexa_Up_Edges(ison,childedge,parentedge)
      use typre
      implicit none
      integer(ip) :: ielty,ison, childedge,parentedge

      if (ison == 1) then
         select case (childedge)
            case (1)
               parentedge = 1
            case(2)
               parentedge = 2
            case(3)
               parentedge = 3
            case default
               parentedge = -1
         end select
      elseif (ison == 2) then
         select case(childedge)
            case(3)
               parentedge = 3
            case(9)
               parentedge = 9
            case(10)
               parentedge = 10
            case default
               parentedge = -1
         end select
      elseif (ison == 3) then
         select case(childedge)
            case(2)
               parentedge = 2
            case(6)
               parentedge = 6
            case(8)
               parentedge = 8
            case default
               parentedge = -1
         end select
       elseif (ison == 4) then
         select case(childedge)
            case(10)
               parentedge = 10
            case(8)
               parentedge = 8
            case(12)
               parentedge = 12
            case default
               parentedge = -1
         end select
      elseif (ison == 5) then
         select case(childedge)
            case(1)
               parentedge = 1
            case(4)
               parentedge = 4
            case(5)
               parentedge = 5
            case default
               parentedge = -1
         end select
      elseif (ison == 6) then
         select case(childedge)
            case(9)
               parentedge = 9
            case(5)
               parentedge = 5
            case(11)
               parentedge = 11
            case default
               parentedge = -1
         end select
      elseif (ison == 7) then
         select case(childedge)
            case(4)
               parentedge = 4
            case(7)
               parentedge = 7
            case(6)
               parentedge = 6
            case default
               parentedge = -1
         end select
      elseif (ison == 8) then
         select case(childedge)
            case(7)
               parentedge = 7
            case(11)
               parentedge = 11
            case(12)
               parentedge = 12
            case default
               parentedge = -1
         end select
      endif
   end subroutine

   subroutine TrilinearHexa_GetNumberofInTheFaceNodesToAdd(nInTheFaceNodes)
      use typre
      implicit none
      integer(ip) :: nInTheFaceNodes

      nInTheFaceNodes = 6

   end subroutine

   subroutine TrilinearHexa_GetFaceForInTheFaceNode(inode,iface)
      use typre
      implicit none
      integer(ip) :: ielty,inode,iface

      iface = inode
   end subroutine

   subroutine TrilinearHexa_GetInTheFacePointFromParent(iface,kchild,knode)
      use typre
      implicit none
      integer(ip) :: iface,kchild,knode

      select case (iface)
         case(1)
            kchild = 1
            knode = 8
         case(2)
            kchild = 5
            knode = 7
         case(3)
            kchild = 1
            knode = 6
         case(4)
            kchild = 3
            knode = 7
         case(5)
            kchild = 1
            knode = 3
         case(6)
            kchild = 2
            knode = 7
      end select

   end subroutine

   subroutine TrilinearHexa_AddNodeToChildInTheFace(ivariation,inode,kpoin,PointerToChildren,pnods,lnods)
      use typre
      implicit none
      integer(ip) :: ivariation,inode,kpoin,PointerToChildren(*), pnods(*),lnods(*)

      integer(ip) :: ielem,ispos

      select case (inode)
         case(1)
            ielem = PointerToChildren(1)
            ispos = pnods(ielem)-1
            lnods(ispos+8) = kpoin

            ielem = PointerToChildren(2)
            ispos = pnods(ielem)-1
            lnods(ispos+4) = kpoin

            ielem = PointerToChildren(3)
            ispos = pnods(ielem)-1
            lnods(ispos+5) = kpoin

            ielem = PointerToChildren(4)
            ispos = pnods(ielem)-1
            lnods(ispos+1) = kpoin
         case(2)
            ielem = PointerToChildren(5)
            ispos = pnods(ielem)-1
            lnods(ispos+7) = kpoin

            ielem = PointerToChildren(6)
            ispos = pnods(ielem)-1
            lnods(ispos+3) = kpoin

            ielem = PointerToChildren(7)
            ispos = pnods(ielem)-1
            lnods(ispos+6) = kpoin

            ielem = PointerToChildren(8)
            ispos = pnods(ielem)-1
            lnods(ispos+2) = kpoin
         case(3)
            ielem = PointerToChildren(1)
            ispos = pnods(ielem)-1
            lnods(ispos+6) = kpoin

            ielem = PointerToChildren(2)
            ispos = pnods(ielem)-1
            lnods(ispos+2) = kpoin

            ielem = PointerToChildren(5)
            ispos = pnods(ielem)-1
            lnods(ispos+5) = kpoin

            ielem = PointerToChildren(6)
            ispos = pnods(ielem)-1
            lnods(ispos+1) = kpoin
         case(4)
            ielem = PointerToChildren(3)
            ispos = pnods(ielem)-1
            lnods(ispos+7) = kpoin

            ielem = PointerToChildren(4)
            ispos = pnods(ielem)-1
            lnods(ispos+3) = kpoin

            ielem = PointerToChildren(7)
            ispos = pnods(ielem)-1
            lnods(ispos+8) = kpoin

            ielem = PointerToChildren(8)
            ispos = pnods(ielem)-1
            lnods(ispos+4) = kpoin

         case(5)
            ielem = PointerToChildren(1)
            ispos = pnods(ielem)-1
            lnods(ispos+3) = kpoin

            ielem = PointerToChildren(3)
            ispos = pnods(ielem)-1
            lnods(ispos+2) = kpoin

            ielem = PointerToChildren(5)
            ispos = pnods(ielem)-1
            lnods(ispos+4) = kpoin

            ielem = PointerToChildren(7)
            ispos = pnods(ielem)-1
            lnods(ispos+1) = kpoin

         case(6)
            ielem = PointerToChildren(2)
            ispos = pnods(ielem)-1
            lnods(ispos+7) = kpoin

            ielem = PointerToChildren(4)
            ispos = pnods(ielem)-1
            lnods(ispos+6) = kpoin

            ielem = PointerToChildren(6)
            ispos = pnods(ielem)-1
            lnods(ispos+8) = kpoin

            ielem = PointerToChildren(8)
            ispos = pnods(ielem)-1
            lnods(ispos+5) = kpoin

      end select
   end subroutine

   subroutine TrilinearHexa_ParentToChildrenElement(ndofn,OldNodalArray,NewNodalArray,ivariation,ichild)
      integer(ip) :: ndofn
      real(rp) :: OldNodalArray(ndofn,8),NewNodalArray(ndofn,8)
      integer(ip) :: ielty,ivariation,ichild

      if (ichild == 1) then
         NewNodalArray(:,1) = OldNodalArray(:,1)
         NewNodalArray(:,2) = (OldNodalArray(:,1)+OldNodalArray(:,2))*0.5
         NewNodalArray(:,3) = (OldNodalArray(:,1)+OldNodalArray(:,2)+OldNodalArray(:,3)+OldNodalArray(:,4))*0.25
         NewNodalArray(:,4) = (OldNodalArray(:,1)+OldNodalArray(:,4))*0.5
         NewNodalArray(:,5) = (OldNodalArray(:,1)+OldNodalArray(:,5))*0.5
         NewNodalArray(:,6) = (OldNodalArray(:,1)+OldNodalArray(:,2)+OldNodalArray(:,5)+OldNodalArray(:,6))*0.25
         NewNodalArray(:,7) = (OldNodalArray(:,1)+OldNodalArray(:,2)+OldNodalArray(:,3)+OldNodalArray(:,4)+OldNodalArray(:,5)+OldNodalArray(:,6)+OldNodalArray(:,7)+OldNodalArray(:,8))*0.125
         NewNodalArray(:,8) = (OldNodalArray(:,1)+OldNodalArray(:,4)+OldNodalArray(:,5)+OldNodalArray(:,8))*0.25

      elseif (ichild == 2) then
         NewNodalArray(:,1) =  (OldNodalArray(:,1)+OldNodalArray(:,5))*0.5
         NewNodalArray(:,2) = (OldNodalArray(:,1)+OldNodalArray(:,2)+OldNodalArray(:,5)+OldNodalArray(:,6))*0.25
         NewNodalArray(:,3) = (OldNodalArray(:,1)+OldNodalArray(:,2)+OldNodalArray(:,3)+OldNodalArray(:,4)+OldNodalArray(:,5)+OldNodalArray(:,6)+OldNodalArray(:,7)+OldNodalArray(:,8))*0.125
         NewNodalArray(:,4) = (OldNodalArray(:,1)+OldNodalArray(:,4)+OldNodalArray(:,5)+OldNodalArray(:,8))*0.25
         NewNodalArray(:,5) = OldNodalArray(:,5)
         NewNodalArray(:,6) =  (OldNodalArray(:,5)+OldNodalArray(:,6))*0.5
         NewNodalArray(:,7) = (OldNodalArray(:,5)+OldNodalArray(:,6)+OldNodalArray(:,7)+OldNodalArray(:,8))*0.25
         NewNodalArray(:,8) =  (OldNodalArray(:,5)+OldNodalArray(:,8))*0.5

      elseif (ichild == 3) then
         NewNodalArray(:,1) = (OldNodalArray(:,1)+OldNodalArray(:,4))*0.5
         NewNodalArray(:,2) = (OldNodalArray(:,1)+OldNodalArray(:,2)+OldNodalArray(:,3)+OldNodalArray(:,4))*0.25
         NewNodalArray(:,3) = (OldNodalArray(:,3)+OldNodalArray(:,4))*0.5
         NewNodalArray(:,4) = OldNodalArray(:,4)
         NewNodalArray(:,5) = (OldNodalArray(:,1)+OldNodalArray(:,4)+OldNodalArray(:,5)+OldNodalArray(:,8))*0.25
         NewNodalArray(:,6) =  (OldNodalArray(:,1)+OldNodalArray(:,2)+OldNodalArray(:,3)+OldNodalArray(:,4)+OldNodalArray(:,5)+OldNodalArray(:,6)+OldNodalArray(:,7)+OldNodalArray(:,8))*0.125
         NewNodalArray(:,7) = (OldNodalArray(:,3)+OldNodalArray(:,4)+OldNodalArray(:,7)+OldNodalArray(:,8))*0.25
         NewNodalArray(:,8) =  (OldNodalArray(:,4)+OldNodalArray(:,8))*0.5
      elseif (ichild == 4) then
         NewNodalArray(:,1) = (OldNodalArray(:,1)+OldNodalArray(:,4)+OldNodalArray(:,5)+OldNodalArray(:,8))*0.25
         NewNodalArray(:,2) = (OldNodalArray(:,1)+OldNodalArray(:,2)+OldNodalArray(:,3)+OldNodalArray(:,4)+OldNodalArray(:,5)+OldNodalArray(:,6)+OldNodalArray(:,7)+OldNodalArray(:,8))*0.125
         NewNodalArray(:,3) = (OldNodalArray(:,3)+OldNodalArray(:,4)+OldNodalArray(:,7)+OldNodalArray(:,8))*0.25
         NewNodalArray(:,4) = (OldNodalArray(:,4)+OldNodalArray(:,8))*0.5
         NewNodalArray(:,5) = (OldNodalArray(:,5)+OldNodalArray(:,8))*0.5
         NewNodalArray(:,6) = (OldNodalArray(:,5)+OldNodalArray(:,6)+OldNodalArray(:,7)+OldNodalArray(:,8))*0.25
         NewNodalArray(:,7) = (OldNodalArray(:,7)+OldNodalArray(:,8))*0.5
         NewNodalArray(:,8) = OldNodalArray(:,8)
      elseif (ichild == 5) then
         NewNodalArray(:,1) = (OldNodalArray(:,1)+OldNodalArray(:,2))*0.5
         NewNodalArray(:,2) = OldNodalArray(:,2)
         NewNodalArray(:,3) = (OldNodalArray(:,2)+OldNodalArray(:,3))*0.5
         NewNodalArray(:,4) = (OldNodalArray(:,1)+OldNodalArray(:,2)+OldNodalArray(:,3)+OldNodalArray(:,4))*0.25
         NewNodalArray(:,5) = (OldNodalArray(:,1)+OldNodalArray(:,2)+OldNodalArray(:,5)+OldNodalArray(:,6))*0.25
         NewNodalArray(:,6) = (OldNodalArray(:,2)+OldNodalArray(:,6))*0.5
         NewNodalArray(:,7) = (OldNodalArray(:,2)+OldNodalArray(:,3)+OldNodalArray(:,6)+OldNodalArray(:,7))*0.25
         NewNodalArray(:,8) = (OldNodalArray(:,1)+OldNodalArray(:,2)+OldNodalArray(:,3)+OldNodalArray(:,4)+OldNodalArray(:,5)+OldNodalArray(:,6)+OldNodalArray(:,7)+OldNodalArray(:,8))*0.125

      elseif (ichild == 6) then
         NewNodalArray(:,1) = (OldNodalArray(:,1)+OldNodalArray(:,2)+OldNodalArray(:,5)+OldNodalArray(:,6))*0.25
         NewNodalArray(:,2) = (OldNodalArray(:,2)+OldNodalArray(:,6))*0.5
         NewNodalArray(:,3) = (OldNodalArray(:,2)+OldNodalArray(:,3)+OldNodalArray(:,6)+OldNodalArray(:,7))*0.25
         NewNodalArray(:,4) = (OldNodalArray(:,1)+OldNodalArray(:,2)+OldNodalArray(:,3)+OldNodalArray(:,4)+OldNodalArray(:,5)+OldNodalArray(:,6)+OldNodalArray(:,7)+OldNodalArray(:,8))*0.125
         NewNodalArray(:,5) = (OldNodalArray(:,5)+OldNodalArray(:,6))*0.5
         NewNodalArray(:,6) =  OldNodalArray(:,6)
         NewNodalArray(:,7) = (OldNodalArray(:,6)+OldNodalArray(:,7))*0.5
         NewNodalArray(:,8) = (OldNodalArray(:,5)+OldNodalArray(:,6)+OldNodalArray(:,7)+OldNodalArray(:,8))*0.25

      elseif (ichild == 7) then
         NewNodalArray(:,1) = (OldNodalArray(:,1)+OldNodalArray(:,2)+OldNodalArray(:,3)+OldNodalArray(:,4))*0.25
         NewNodalArray(:,2) = (OldNodalArray(:,2)+OldNodalArray(:,3))*0.5
         NewNodalArray(:,3) = OldNodalArray(:,3)
         NewNodalArray(:,4) = (OldNodalArray(:,3)+OldNodalArray(:,4))*0.5
         NewNodalArray(:,5) = (OldNodalArray(:,1)+OldNodalArray(:,2)+OldNodalArray(:,3)+OldNodalArray(:,4)+OldNodalArray(:,5)+OldNodalArray(:,6)+OldNodalArray(:,7)+OldNodalArray(:,8))*0.125
         NewNodalArray(:,6) = (OldNodalArray(:,2)+OldNodalArray(:,3)+OldNodalArray(:,6)+OldNodalArray(:,7))*0.25
         NewNodalArray(:,7) = (OldNodalArray(:,3)+OldNodalArray(:,7))*0.5
         NewNodalArray(:,8) = (OldNodalArray(:,3)+OldNodalArray(:,4)+OldNodalArray(:,7)+OldNodalArray(:,8))*0.25
      elseif (ichild == 8) then
         NewNodalArray(:,1) = (OldNodalArray(:,1)+OldNodalArray(:,2)+OldNodalArray(:,3)+OldNodalArray(:,4)+OldNodalArray(:,5)+OldNodalArray(:,6)+OldNodalArray(:,7)+OldNodalArray(:,8))*0.125
         NewNodalArray(:,2) = (OldNodalArray(:,2)+OldNodalArray(:,3)+OldNodalArray(:,6)+OldNodalArray(:,7))*0.25
         NewNodalArray(:,3) = (OldNodalArray(:,3)+OldNodalArray(:,7))*0.5
         NewNodalArray(:,4) = (OldNodalArray(:,3)+OldNodalArray(:,4)+OldNodalArray(:,7)+OldNodalArray(:,8))*0.25
         NewNodalArray(:,5) = (OldNodalArray(:,5)+OldNodalArray(:,6)+OldNodalArray(:,7)+OldNodalArray(:,8))*0.25
         NewNodalArray(:,6) = (OldNodalArray(:,6)+OldNodalArray(:,7))*0.5
         NewNodalArray(:,7) = OldNodalArray(:,7)
         NewNodalArray(:,8) = (OldNodalArray(:,7)+OldNodalArray(:,8))*0.5
      endif

   end subroutine

   subroutine TrilinearHexa_ChildrenToParentElement(ndofn,OldNodalArray,NewNodalArray,ivariation,ichild)
      integer(ip) :: ndofn
      real(rp) :: OldNodalArray(ndofn,8),NewNodalArray(ndofn,8)
      integer(ip) :: ivariation,ichild

      if (ichild == 1) then
         NewNodalArray(:,1) = OldNodalArray(:,1)
      elseif (ichild == 2) then
         NewNodalArray(:,5) = OldNodalArray(:,5)
      elseif (ichild == 3) then
         NewNodalArray(:,4) = OldNodalArray(:,4)
      elseif (ichild == 4) then
         NewNodalArray(:,8) = OldNodalArray(:,8)
      elseif (ichild == 5) then
         NewNodalArray(:,2) = OldNodalArray(:,2)
      elseif (ichild == 6) then
         NewNodalArray(:,6) = OldNodalArray(:,6)
      elseif (ichild == 7) then
         NewNodalArray(:,3) = OldNodalArray(:,3)
      elseif (ichild == 8) then
         NewNodalArray(:,7) = OldNodalArray(:,7)
      endif

   end subroutine

end module

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


! MD5 of template: 888972e9e180574c8c83291b6532a6de
! Random number generation related routines
!
! ------------------------------------------------------------------------------
! Copyright (c) 2009-13, Thomas P. Robitaille
!
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
!  * Redistributions of source code must retain the above copyright notice, this
!    list of conditions and the following disclaimer.
!
!  * Redistributions in binary form must reproduce the above copyright notice,
!    this list of conditions and the following disclaimer in the documentation
!    and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
! FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ------------------------------------------------------------------------------

module lib_random

  implicit none
  save

  private

  integer,parameter :: sp = selected_real_kind(p=6,r=37)
  integer,parameter :: dp = selected_real_kind(p=15,r=307)

  integer :: idum = -1204132124
  real(dp) :: u(97)

  real(dp),parameter :: pi = 3.14159265358979323846_dp
  real(sp),parameter :: pi_sp = 3.14159265358979323846_sp
  real(dp),parameter :: pi_dp = 3.14159265358979323846_dp

  real(dp),parameter :: twopi = 2._dp * pi
  real(sp),parameter :: twopi_sp = pi_sp + pi_sp
  real(dp),parameter :: twopi_dp = pi_dp + pi_dp

  public :: set_seed
  public :: set_seed_64

  public :: random
  interface random
     module procedure random_sp
     module procedure random_dp
  end interface random

  public :: random_exp
  interface random_exp
     module procedure random_exp_sp
     module procedure random_exp_dp
  end interface random_exp

  public :: random_uni
  interface random_uni
     module procedure random_uni_sp
     module procedure random_uni_dp
  end interface random_uni

  public :: random_gau
  interface random_gau
     module procedure random_gau_sp
     module procedure random_gau_dp
  end interface random_gau

  public :: random_sphere
  interface random_sphere
     module procedure random_sphere_sp
     module procedure random_sphere_dp
  end interface random_sphere

  public :: random_poisson
  interface random_poisson
     module procedure random_poisson_sp
     module procedure random_poisson_dp
  end interface random_poisson

  public :: random_planck_frequency
  interface random_planck_frequency
     module procedure random_planck_frequency_sp
     module procedure random_planck_frequency_dp
  end interface random_planck_frequency

contains

  subroutine set_seed(seed)
    implicit none
    integer,intent(in) :: seed
    idum = -abs(seed)
    call set_seed_64(abs(seed), 987654321)
  end subroutine set_seed

  subroutine set_seed_64(seed1,seed2)
    implicit none
    integer,intent(in) :: seed1,seed2
    integer i,j,x,y
    real(dp) :: s,t
    x=seed1
    y=seed2
    do i=1,97
       s=0._dp
       t=0.5_dp
       do j=1,53
          x=mod(6969*x,65543)
          y=mod(8888*x,65579)
          if (iand(ieor(x,y),32).gt.0) s=s+t
          t=0.5_dp*t
       end do
       u(i)=s
    end do
  end subroutine set_seed_64

  subroutine random_string(string)
    implicit none
    character(len=*),intent(out) :: string
    integer :: i,j
    real(dp) :: xi
    do i=1,len(string)
       call random(xi)
       j = int(xi*26._dp)
       call random(xi)
       if(xi < 0.5) then
          j = j + 65
       else
          j = j + 97
       end if
       string(i:i) = char(j)
    end do
  end subroutine random_string

  subroutine random_sp(xi)
    ! Random number between 0 and 1
    ! Inspired by Numerical Recipes
    implicit none
    real(sp),intent(out) :: xi
    real(sp),save        :: am
    integer, parameter :: ia=16807,im=2147483647,iq=127773,ir=2836
    integer, save :: ix=-1,iy=-1,k
    if (idum <= 0 .or. iy < 0) then
       am=nearest(1.0_sp,-1.0_sp)/im
       iy=ior(ieor(888889999,abs(idum)),1)
       ix=ieor(777755555,abs(idum))
       idum=abs(idum)+1
    end if
    ix=ieor(ix,ishft(ix,13))
    ix=ieor(ix,ishft(ix,-17))
    ix=ieor(ix,ishft(ix,5))
    k=iy/iq
    iy=ia*(iy-k*iq)-ir*k
    if (iy < 0) iy=iy+im
    xi=am*ior(iand(im,ieor(ix,iy)),1)
  end subroutine random_sp

  subroutine random_dp(xi)
    ! Random number between 0 and 1
    ! Based on "The 64-bit universal RNG", Marsaglia & Tsang (2004)
    implicit none
    real(dp),intent(out) :: xi
    integer,save :: i=97
    integer,save :: j=33
    real(dp),save :: c=0
    real(dp) :: x
    real(dp), parameter :: r=9007199254740881._dp/9007199254740992._dp
    real(dp), parameter :: d=362436069876._dp/9007199254740992._dp
    x=u(i)-u(j)
    if (x.lt.0.0) x=x+1.0_dp
    u(i)=x
    i=i-1
    if (i.eq.0) i=97
    j=j-1
    if (j.eq.0) j=97
    c=c-d
    if (c.lt.0.0) c=c+r
    x=x-c
    xi=x
    if (x.lt.0.) xi=x+1._dp
    return
  end subroutine random_dp

  subroutine random_uni_dp(xi,a,b)
    ! Uniform random number between a and b
    implicit none
    real(dp),intent(in)  :: a,b
    real(dp),intent(out) :: xi
    call random(xi)
    xi = a + (b-a) * xi
  end subroutine random_uni_dp

  subroutine random_gau_dp(xi,c,w)
    ! Gaussian random number with center c and 1-sigma w
    implicit none
    real(dp),intent(in)  :: c,w
    real(dp),intent(out) :: xi
    real(dp) :: r,g1,g2,x,y,t
    do
       call random_uni(x,-1._dp,+1._dp)
       call random_uni(y,-1._dp,+1._dp)
       r = x*x + y*y
       if(r.lt.1) exit
    end do
    t=sqrt(-2._dp*log(r)/r)
    g1 = x*t
    g2 = y*t
    xi = g1*w+c
  end subroutine random_gau_dp

  subroutine random_exp_dp(xi)
    ! Random number sampled from exp(-tau)
    implicit none
    real(dp),intent(out) :: xi
    do
       call random(xi)
       if(xi < 1._dp) exit
    end do
    xi = - log( 1._dp - xi )
  end subroutine random_exp_dp

  subroutine random_sphere_dp(mu,phi)
    ! Random longitude/latitude on a sphere
    implicit none
    real(dp),intent(out) :: mu,phi
    call random_uni(mu,-1._dp,+1._dp)
    call random_uni(phi,0._dp,twopi_dp)
  end subroutine random_sphere_dp

  subroutine random_poisson_dp(xi,xm)
    ! Random number drawn from Poisson distribution with
    ! expected value xm

    implicit none

    real(dp),intent(in) :: xm
    real(dp),intent(out) :: xi
    real(dp) :: em,harvest,t,y

    real(dp), save :: alxm,g,oldm=-1.0_dp,sq

    if (xm < 12.0) then
       if (xm /= oldm) then
          oldm=xm
          g=exp(-xm)
       end if
       em=-1
       t=1.0
       do
          em=em+1.0_dp
          call random_dp(harvest)
          t=t*harvest
          if (t <= g) exit
       end do
    else
       if (xm /= oldm) then
          oldm=xm
          sq=sqrt(2.0_dp*xm)
          alxm=log(xm)
          g=xm*alxm-gammln_dp(xm+1.0_dp)
       end if
       do
          do
             call random_dp(harvest)
             y=tan(pi*harvest)
             em=sq*y+xm
             if (em >= 0.0) exit
          end do
          em=int(em)
          t=0.9_dp*(1.0_dp+y**2)*exp(em*alxm-gammln_dp(em+1.0_dp)-g)
          call random_dp(harvest)
          if (harvest <= t) exit
       end do
    end if
    xi = em
  end subroutine random_poisson_dp

  subroutine random_planck_frequency_dp(nu,T)

    ! Random frequency sampled from a planck function with temperature T

    ! The algorithm is taken from 'Sampling a random variable distributed
    ! according to planck's law' by Barnett and Canfield

    implicit none

    real(dp),intent(in) :: T
    real(dp),intent(out) :: nu
    real(dp) :: x,r,r1,r2,r3,r4,a,y,z

    real(dp),parameter :: k  = 1.3806503e-23_dp ! J/K
    real(dp),parameter :: h  = 6.626068e-34_dp ! J.s

    ! Sample a random number from x^3/(exp(x)-1)

    do

       call random(r1)
       call random(r2)
       call random(r3)
       call random(r4)

       r = r1*r2*r3*r4

       if(r > 0._dp) exit

    end do

    x = - log(r)

    a = 1._dp
    y = 1._dp
    z = 1._dp

    call random(r1)
    do
       if(1.08232_dp*r1 <= a) exit
       y = y + 1._dp
       z = 1._dp/y
       a = a + z*z*z*z
    end do
    x = x * z

    ! Convert to frequency

    nu = x * k * T / h

  end subroutine random_planck_frequency_dp

  real(dp) function gammln_dp(xx)

    implicit none

    real(dp),intent(in) :: xx
    integer :: j
    real(dp) :: ser,tmp,x,y

    real(dp),save :: cof(6) = (/76.18009172947146_dp,&
         &-86.50532032941677_dp,24.01409824083091_dp,&
         &-1.231739572450155_dp,.1208650973866179e-2_dp,&
         &-.5395239384953e-5_dp/)

    real(dp),save :: stp = 2.5066282746310005_dp

    x=xx
    y=x
    tmp=x+5.5_dp
    tmp=(x+0.5_dp)*log(tmp)-tmp
    ser=1.000000000190015_dp

    do j=1,6
       y=y+1._dp
       ser=ser+cof(j)/y
    end do

    gammln_dp=tmp+log(stp*ser/x)
    return

  end function gammln_dp

  subroutine random_uni_sp(xi,a,b)
    ! Uniform random number between a and b
    implicit none
    real(sp),intent(in)  :: a,b
    real(sp),intent(out) :: xi
    call random(xi)
    xi = a + (b-a) * xi
  end subroutine random_uni_sp

  subroutine random_gau_sp(xi,c,w)
    ! Gaussian random number with center c and 1-sigma w
    implicit none
    real(sp),intent(in)  :: c,w
    real(sp),intent(out) :: xi
    real(sp) :: r,g1,g2,x,y,t
    do
       call random_uni(x,-1._sp,+1._sp)
       call random_uni(y,-1._sp,+1._sp)
       r = x*x + y*y
       if(r.lt.1) exit
    end do
    t=sqrt(-2._sp*log(r)/r)
    g1 = x*t
    g2 = y*t
    xi = g1*w+c
  end subroutine random_gau_sp

  subroutine random_exp_sp(xi)
    ! Random number sampled from exp(-tau)
    implicit none
    real(sp),intent(out) :: xi
    do
       call random(xi)
       if(xi < 1._sp) exit
    end do
    xi = - log( 1._sp - xi )
  end subroutine random_exp_sp

  subroutine random_sphere_sp(mu,phi)
    ! Random longitude/latitude on a sphere
    implicit none
    real(sp),intent(out) :: mu,phi
    call random_uni(mu,-1._sp,+1._sp)
    call random_uni(phi,0._sp,twopi_sp)
  end subroutine random_sphere_sp

  subroutine random_poisson_sp(xi,xm)
    ! Random number drawn from Poisson distribution with
    ! expected value xm

    implicit none

    real(sp),intent(in) :: xm
    real(sp),intent(out) :: xi
    real(sp) :: em,harvest,t,y

    real(sp), save :: alxm,g,oldm=-1.0_sp,sq

    if (xm < 12.0) then
       if (xm /= oldm) then
          oldm=xm
          g=exp(-xm)
       end if
       em=-1
       t=1.0
       do
          em=em+1.0_sp
          call random_sp(harvest)
          t=t*harvest
          if (t <= g) exit
       end do
    else
       if (xm /= oldm) then
          oldm=xm
          sq=sqrt(2.0_sp*xm)
          alxm=log(xm)
          g=xm*alxm-gammln_sp(xm+1.0_sp)
       end if
       do
          do
             call random_sp(harvest)
             y=tan(pi*harvest)
             em=sq*y+xm
             if (em >= 0.0) exit
          end do
          em=int(em)
          t=0.9_sp*(1.0_sp+y**2)*exp(em*alxm-gammln_sp(em+1.0_sp)-g)
          call random_sp(harvest)
          if (harvest <= t) exit
       end do
    end if
    xi = em
  end subroutine random_poisson_sp

  subroutine random_planck_frequency_sp(nu,T)

    ! Random frequency sampled from a planck function with temperature T

    ! The algorithm is taken from 'Sampling a random variable distributed
    ! according to planck's law' by Barnett and Canfield

    implicit none

    real(sp),intent(in) :: T
    real(sp),intent(out) :: nu
    real(sp) :: x,r,r1,r2,r3,r4,a,y,z

    real(sp),parameter :: k  = 1.3806503e-23_sp ! J/K
    real(sp),parameter :: h  = 6.626068e-34_sp ! J.s

    ! Sample a random number from x^3/(exp(x)-1)

    do

       call random(r1)
       call random(r2)
       call random(r3)
       call random(r4)

       r = r1*r2*r3*r4

       if(r > 0._dp) exit

    end do

    x = - log(r)

    a = 1._sp
    y = 1._sp
    z = 1._sp

    call random(r1)
    do
       if(1.08232_sp*r1 <= a) exit
       y = y + 1._sp
       z = 1._sp/y
       a = a + z*z*z*z
    end do
    x = x * z

    ! Convert to frequency

    nu = x * k * T / h

  end subroutine random_planck_frequency_sp

  real(sp) function gammln_sp(xx)

    implicit none

    real(sp),intent(in) :: xx
    integer :: j
    real(sp) :: ser,tmp,x,y

    real(sp),save :: cof(6) = (/76.18009172947146_sp,&
         &-86.50532032941677_sp,24.01409824083091_sp,&
         &-1.231739572450155_sp,.1208650973866179e-2_sp,&
         &-.5395239384953e-5_sp/)

    real(sp),save :: stp = 2.5066282746310005_sp

    x=xx
    y=x
    tmp=x+5.5_sp
    tmp=(x+0.5_sp)*log(tmp)-tmp
    ser=1.000000000190015_sp

    do j=1,6
       y=y+1._sp
       ser=ser+cof(j)/y
    end do

    gammln_sp=tmp+log(stp*ser/x)
    return

  end function gammln_sp

end module lib_random


text:
C
C NLP-UNCONSTRAINED-CLI/NELDER-MEAD/__ORIG/SRC/NELMIN.F
C =============================================================================
C NONLINEAR OPTIMIZATION ALGORITHMS MULTILANG. VERSION 0.1
C =============================================================================
C NONLINEAR PROGRAMMING ALGORITHMS AS THE (UN-)CONSTRAINED MINIMIZATION
C PROBLEMS WITH THE FOCUS ON THEIR NUMERICAL EXPRESSION USING VARIOUS
C PROGRAMMING LANGUAGES.
C
C THIS IS THE NELDER-MEAD NONLINEAR UNCONSTRAINED MINIMIZATION ALGORITHM.
C =============================================================================
C WRITTEN BY RADISLAV (RADICCHIO) GOLUBTSOV, 2015-2020
C
C THIS IS FREE AND UNENCUMBERED SOFTWARE RELEASED INTO THE PUBLIC DOMAIN.
C
C ANYONE IS FREE TO COPY, MODIFY, PUBLISH, USE, COMPILE, SELL, OR
C DISTRIBUTE THIS SOFTWARE, EITHER IN SOURCE CODE FORM OR AS A COMPILED
C BINARY, FOR ANY PURPOSE, COMMERCIAL OR NON-COMMERCIAL, AND BY ANY
C MEANS.
C
C (SEE THE LICENSE FILE AT THE TOP OF THE SOURCE TREE.)

C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C *                                                                           *
C *             SEE THE NELMIN.TXT.C FILE FOR DETAILS DESCRIPTION.            *
C *                                                                           *
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

C     === THE USER-SUPPLIED OBJECTIVE FUNCTION F(X).
      FUNCTION F(X)
          IMPLICIT NONE

C         CONSTANT. THE MAXIMUM NUMBER OF VARIABLES.
          INTEGER    VARS
          PARAMETER (VARS = 20)

C         RETURN VAL. THE OBJECTIVE FUNCTION VALUE.
          DOUBLE PRECISION F

C         ARG. THE POINT AT WHICH F(X) SHOULD BE EVALUATED.
          DOUBLE PRECISION X(VARS)

#ifndef WOODS
C     ROSENBROCK'S CLASSIC PARABOLIC VALLEY ('BANANA') FUNCTION.
          DOUBLE PRECISION A
          DOUBLE PRECISION B

          A = X(2) - X(1) ** 2
          B = 1    - X(1)

          F = 100 * (A ** 2) + (B ** 2)
#else
C     WOODS -- A LA MORE, GARBOW AND HILLSTROM (TOMS ALGORITHM 566).
          DOUBLE PRECISION S1
          DOUBLE PRECISION S2
          DOUBLE PRECISION S3
          DOUBLE PRECISION T1
          DOUBLE PRECISION T2
          DOUBLE PRECISION T3
          DOUBLE PRECISION T4
          DOUBLE PRECISION T5

          S1 = X(2) - X(1) ** 2
          S2 = 1    - X(1)
          S3 = X(2) - 1

          T1 = X(4) - X(3) ** 2
          T2 = 1    - X(3)
          T3 = X(4) - 1

          T4 = S3 + T3
          T5 = S3 - T3

          F = 100 * (S1 ** 2) + (S2 ** 2)
     *       + 90 * (T1 ** 2) + (T2 ** 2)
     *       + 10 * (T4 ** 2) + (T5 ** 2) / 10
#endif
          RETURN
      END

C     === MAIN OPTIMIZATION SUBROUTINE.
C     THE NELMIN SUBROUTINE ITSELF (NELDER-MEAD MINIMIZATION).
      SUBROUTINE NELMIN(   F,      N,  START,   XMIN, YNEWLO, REQMIN,
     *                  STEP, KONVGE, KCOUNT, ICOUNT, NUMRES, IFAULT)

          IMPLICIT NONE

C         CONSTANT. THE MAXIMUM NUMBER OF VARIABLES.
          INTEGER    VARS
          PARAMETER (VARS = 20)

C         CONSTANT. THE REFLECTION COEFFICIENT.
          DOUBLE PRECISION RCOEFF
          PARAMETER       (RCOEFF = 1.0D+00)

C         CONSTANT. THE EXTENSION COEFFICIENT.
          DOUBLE PRECISION ECOEFF
          PARAMETER       (ECOEFF = 2.0D+00)

C         CONSTANT. THE CONTRACTION COEFFICIENT.
          DOUBLE PRECISION CCOEFF
          PARAMETER       (CCOEFF = 0.5D+00)

C         CONSTANT. THE OPTIMALITY FACTOR.
          DOUBLE PRECISION EPS
          PARAMETER       (EPS = 0.001D+00)

C         ARG. THE OBJECTIVE FUNCTION F(X).
          DOUBLE PRECISION F
          EXTERNAL         F

C         ARG. THE NUMBER OF VARIABLES.
          INTEGER N

C         ARG. THE STARTING POINT FOR THE ITERATION.
          DOUBLE PRECISION START(VARS)

C         ARG. THE COORDINATES OF THE POINT WHICH IS ESTIMATED
C              TO MINIMIZE THE FUNCTION.
          DOUBLE PRECISION XMIN(VARS)

C         ARG. THE MINIMUM VALUE OF THE FUNCTION.
          DOUBLE PRECISION YNEWLO

C         ARG. THE TERMINATING LIMIT FOR THE VARIANCE
C              OF FUNCTION VALUES.
          DOUBLE PRECISION REQMIN

C         ARG. THE SIZE AND SHAPE OF THE INITIAL SIMPLEX.
          DOUBLE PRECISION STEP(VARS)

C         ARG. THE CONVERGENCE CHECK.
          INTEGER KONVGE

C         ARG. THE MAXIMUM NUMBER OF FUNCTION EVALUATIONS.
          INTEGER KCOUNT

C         ARG. THE NUMBER OF FUNCTION EVALUATIONS USED.
          INTEGER ICOUNT

C         ARG. THE NUMBER OF RESTARTS.
          INTEGER NUMRES

C         ARG. THE ERROR INDICATOR.
          INTEGER IFAULT

          INTEGER JCOUNT
          INTEGER NN
          INTEGER I
          INTEGER J
          INTEGER ILO
          INTEGER IHI
          INTEGER L

          DOUBLE PRECISION DN
          DOUBLE PRECISION DNN
          DOUBLE PRECISION DEL
          DOUBLE PRECISION RQ
          DOUBLE PRECISION P(VARS, VARS + 1)
          DOUBLE PRECISION Y(VARS + 1)
          DOUBLE PRECISION X
          DOUBLE PRECISION YLO
          DOUBLE PRECISION Z
          DOUBLE PRECISION PBAR(VARS)
          DOUBLE PRECISION PSTAR(VARS)
          DOUBLE PRECISION YSTAR
          DOUBLE PRECISION P2STAR(VARS)
          DOUBLE PRECISION Y2STAR

C         CHECK THE INPUT PARAMETERS.
          IF (REQMIN .LE. 0.0D+00) THEN
              IFAULT = 1

              RETURN
          END IF

          IF (N .LT. 1) THEN
              IFAULT = 1

              RETURN
          END IF

          IF (VARS .LT. N) THEN
              IFAULT = 1

              RETURN
          END IF

          IF (KONVGE .LT. 1) THEN
              IFAULT = 1

              RETURN
          END IF

          ICOUNT = 0
          NUMRES = 0

          JCOUNT = KONVGE
          DN     = DBLE(N)
          NN     = N + 1
          DNN    = DBLE(NN)
          DEL    = 1.0D+00
          RQ     = REQMIN * DN

C         CONSTRUCTION OF INITIAL SIMPLEX.
1000      CONTINUE

          DO 10 I = 1, N
              P(I, NN) = START(I)
10        CONTINUE

          Y(NN) = F(START)

          DO 20 J = 1, N
              X        = START(J)
              START(J) = START(J) + STEP(J) * DEL

              DO 30 I = 1, N
                  P(I, J) = START(I)
30            CONTINUE

              Y(J) = F(START)

              START(J) = X
20        CONTINUE

          ICOUNT = ICOUNT + NN

C         THE SIMPLEX CONSTRUCTION IS COMPLETE.

C         FIND HIGHEST AND LOWEST Y VALUES.
C         YNEWLO = Y(IHI) INDICATES THE VERTEX OF THE SIMPLEX
C         TO BE REPLACED.
          YLO = Y(1)
          ILO =   1

          DO 40 I = 2, NN
              IF (Y(I) .LT. YLO) THEN
                  YLO = Y(I)
                  ILO = I
              END IF
40        CONTINUE

2000      CONTINUE

          YNEWLO = Y(1)
          IHI    =   1

          DO 50 I = 2, NN
              IF (YNEWLO .LT. Y(I)) THEN
                  YNEWLO = Y(I)
                  IHI    = I
              END IF
50        CONTINUE

C         CALCULATE PBAR, THE CENTROID OF THE SIMPLEX VERTICES
C         EXCEPTING THE VERTEX WITH Y VALUE YNEWLO.
          DO 60 I = 1, N
              Z = 0.0D+00

              DO 70 J = 1, NN
                  Z = Z + P(I, J)
70            CONTINUE

              Z       = Z - P(I, IHI)
              PBAR(I) = Z / DN
60        CONTINUE

C         REFLECTION THROUGH THE CENTROID.
          DO 80 I = 1, N
              PSTAR(I) = PBAR(I) + RCOEFF * (PBAR(I) - P(I, IHI))
80        CONTINUE

          YSTAR = F(PSTAR)

          ICOUNT = ICOUNT + 1

C         SUCCESSFUL REFLECTION, SO EXTENSION.
          IF (YSTAR .LT. YLO) THEN
              DO 90 I = 1, N
                  P2STAR(I) = PBAR(I) + ECOEFF * (PSTAR(I) - PBAR(I))
90            CONTINUE

              Y2STAR = F(P2STAR)

              ICOUNT = ICOUNT + 1

C             CHECK EXTENSION.
              IF (YSTAR .LT. Y2STAR) THEN
                  DO 100 I = 1, N
                      P(I, IHI) = PSTAR(I)
100               CONTINUE

                  Y(IHI) = YSTAR
C             RETAIN EXTENSION OR CONTRACTION.
              ELSE
                  DO 110 I = 1, N
                      P(I, IHI) = P2STAR(I)
110               CONTINUE

                  Y(IHI) = Y2STAR
              END IF
C         NO EXTENSION.
          ELSE
              L = 0

              DO 120 I = 1, NN
                  IF (YSTAR .LT. Y(I)) THEN
                      L = L + 1
                  END IF
120           CONTINUE

              IF (1 .LT. L) THEN
                  DO 130 I = 1, N
                      P(I, IHI) = PSTAR(I)
130               CONTINUE

                  Y(IHI) = YSTAR
C             CONTRACTION ON THE Y(IHI) SIDE OF THE CENTROID.
              ELSE IF (L .EQ. 0) THEN
                  DO 140 I = 1, N
                      P2STAR(I) = PBAR(I)
     *                          + CCOEFF * (P(I, IHI) - PBAR(I))
140               CONTINUE

                  Y2STAR = F(P2STAR)

                  ICOUNT = ICOUNT + 1

C                 CONTRACT THE WHOLE SIMPLEX.
                  IF (Y(IHI) .LT. Y2STAR) THEN
                      DO 150 J = 1, NN
                          DO 160 I = 1, N
                              P(I, J) = (P(I, J) + P(I, ILO)) * CCOEFF
                              XMIN(I) =  P(I, J)
160                       CONTINUE

                          Y(J) = F(XMIN)
150                   CONTINUE

                      ICOUNT = ICOUNT + NN

                      IF (KCOUNT .LT. ICOUNT) THEN
                          GO TO 3000
                      END IF

                      YLO = Y(1)
                      ILO =   1

                      DO 170 I = 2, NN
                          IF (Y(I) .LT. YLO) THEN
                              YLO = Y(I)
                              ILO = I
                          END IF
170                   CONTINUE

                      GO TO 2000
C                 RETAIN CONTRACTION.
                  ELSE
                      DO 180 I = 1, N
                          P(I, IHI) = P2STAR(I)
180                   CONTINUE

                      Y(IHI) = Y2STAR
                  END IF
C             CONTRACTION ON THE REFLECTION SIDE OF THE CENTROID.
              ELSE IF (L .EQ. 1) THEN
                  DO 190 I = 1, N
                      P2STAR(I) = PBAR(I)
     *                          + CCOEFF * (PSTAR(I) - PBAR(I))
190               CONTINUE

                  Y2STAR = F(P2STAR)

                  ICOUNT = ICOUNT + 1

C                 RETAIN REFLECTION?
                  IF (Y2STAR .LE. YSTAR) THEN
                      DO 200 I = 1, N
                          P(I, IHI) = P2STAR(I)
200                   CONTINUE

                      Y(IHI) = Y2STAR
                  ELSE
                      DO 210 I = 1, N
                          P(I, IHI) = PSTAR(I)
210                   CONTINUE

                      Y(IHI) = YSTAR
                  END IF
              END IF
          END IF

C         CHECK IF YLO IMPROVED.
          IF (Y(IHI) .LT. YLO) THEN
              YLO = Y(IHI)
              ILO = IHI
          END IF

          JCOUNT = JCOUNT - 1

          IF (JCOUNT .NE. 0) THEN
              GO TO 2000
          END IF

C         CHECK TO SEE IF MINIMUM REACHED.
          IF (ICOUNT .LE. KCOUNT) THEN
              JCOUNT = KONVGE
              Z      = 0.0D+00

              DO 220 I = 1, NN
                  Z = Z + Y(I)
220           CONTINUE

              X = Z / DNN
              Z = 0.0D+00

              DO 230 I = 1, NN
                  Z = Z + (Y(I) - X) ** 2
230           CONTINUE

              IF (RQ .LT. Z) THEN
                  GO TO 2000
              END IF
          END IF

C         FACTORIAL TESTS TO CHECK THAT YNEWLO IS A LOCAL MINIMUM.
3000      CONTINUE

          DO 240 I = 1, N
              XMIN(I) = P(I, ILO)
240       CONTINUE

          YNEWLO = Y(ILO)

          IF (KCOUNT .LT. ICOUNT) THEN
              IFAULT = 2

              RETURN
          END IF

          IFAULT = 0

          DO 250 I = 1, N
              DEL     = STEP(I) * EPS
              XMIN(I) = XMIN(I) + DEL

              Z = F(XMIN)

              ICOUNT = ICOUNT + 1

              IF (Z .LT. YNEWLO) THEN
                  IFAULT = 2

                  GO TO 4000
              END IF

              XMIN(I) = XMIN(I) - (DEL * 2)

              Z = F(XMIN)

              ICOUNT = ICOUNT + 1

              IF (Z .LT. YNEWLO) THEN
                  IFAULT = 2

                  GO TO 4000
              END IF

              XMIN(I) = XMIN(I) + DEL
250       CONTINUE

4000      CONTINUE

          IF (IFAULT .EQ. 0) THEN
              RETURN
          END IF

C         RESTART THE PROCEDURE.
          DO 260 I = 1, N
              START(I) = XMIN(I)
260       CONTINUE

          DEL    = EPS
          NUMRES = NUMRES + 1

          GO TO 1000
      END

C     === MAIN PROGRAM.
      PROGRAM AMOEBA
          IMPLICIT NONE

C         CONSTANT. THE MAXIMUM NUMBER OF VARIABLES.
          INTEGER    VARS
          PARAMETER (VARS = 20)

          INTEGER N
          INTEGER KONVGE
          INTEGER KCOUNT
          INTEGER I
          INTEGER ICOUNT
          INTEGER NUMRES
          INTEGER IFAULT

          DOUBLE PRECISION START(VARS)
          DOUBLE PRECISION REQMIN
          DOUBLE PRECISION STEP(VARS)
          DOUBLE PRECISION YNEWLO
          DOUBLE PRECISION XMIN(VARS)

C         PROTO REF. THE OBJECTIVE FUNCTION F(X).
          DOUBLE PRECISION F

#ifndef WOODS
C     STARTING GUESS FOR ROSENBROCK'S TEST FUNCTION.
          PRINT 10
10        FORMAT (/,   'TEST01',
     *            /2X, 'Apply NELMIN to ROSENBROCK function.')

          N        =  2
          START(1) = -1.2
          START(2) =  1.0
#else
C     STARTING GUESS TEST PROBLEM 'WOODS'.
          PRINT 50
50        FORMAT (/,   'TEST05',
     *            /2X, 'Apply NELMIN to WOODS function.')

          N        =  4
          START(1) = -3.0
          START(2) = -1.0
          START(3) = -3.0
          START(4) = -1.0
#endif

          REQMIN  = 1.0D-08
          STEP(1) = 1.0
          STEP(2) = 1.0

#ifdef WOODS
          STEP(3) = 1.0
          STEP(4) = 1.0
#endif

C         KONVGE AND KCOUNT HAVE SAME VALUES FOR ALL TEST PRBMS:
C         SHOULD THEY BE VARIED?
          KONVGE  = 10
          KCOUNT  = 500

          PRINT 60
60        FORMAT (/2X, 'Starting point X:', /)

          DO 70 I = 1, N
              PRINT 80, START(I), START(I)
80            FORMAT (            2X, G14.6, 12X, '->', 4X, F20.12)
70        CONTINUE

          YNEWLO = F(START)

          PRINT 90, YNEWLO, YNEWLO
90        FORMAT (/2X, 'F(X)  =', 1X, G14.6,  4X, '->', 4X, F20.12)

          CALL NELMIN(   F,      N,  START,   XMIN, YNEWLO, REQMIN,
     *                STEP, KONVGE, KCOUNT, ICOUNT, NUMRES, IFAULT)

          PRINT 100, IFAULT
100       FORMAT (/2X, 'Return code IFAULT   =', 1X, I8)

          PRINT 110
110       FORMAT (/2X, 'Estimate of minimizing value X*:', /)

          DO 120 I = 1, N
              PRINT 130, XMIN(I), XMIN(I)
130           FORMAT (            2X, G14.6, 12X, '->', 4X, F20.12)
120       CONTINUE

          PRINT 140, YNEWLO, YNEWLO
140       FORMAT (/2X, 'F(X*) =', 1X, G14.6,  4X, '->', 4X, F20.12)

          PRINT 150, ICOUNT
150       FORMAT (/2X, 'Number of iterations =', 1X, I8)

          PRINT 160, NUMRES
160       FORMAT ( 2X, 'Number of restarts   =', 1X, I8)
      END


text:
        implicit real *8 (a-h,o-z)
        complex *16 z,j0,ima,ztest(1000 000),jtest(1000 000),
     1    jtest2(1000 000),jtest3(1000 000)
        real *8 rtest(1000 0000),thetatest(1000 0000),jtestr(1000 0000),
     1      jtesti(1000 0000),xtest(1000 0000),ytest(1000 0000),
     2      rerrors(5,4),times(4,2)
        data ima/(0.0d0,1.0d0)/

        call prini(6,13)

        do i=1,4
           do j=1,5
              rerrors(j,i)=0
           enddo
        enddo

        ntest=10000
        OPEN(UNIT=66,FILE="../mathematica/rm1test.txt",
     1      FORM="FORMATTED",STATUS="OLD",ACTION="READ")
        do i=1,ntest
           read(66,*) rtest(i)
        enddo
        close(66)

        OPEN(UNIT=66,FILE="../mathematica/thetam1test.txt",
     1      FORM="FORMATTED",STATUS="OLD",ACTION="READ")
        do i=1,ntest
           read(66,*) thetatest(i)
        enddo
        close(66)

        OPEN(UNIT=66,FILE="../mathematica/jm1testr.txt",
     1      FORM="FORMATTED",STATUS="OLD",ACTION="READ")
        do i=1,ntest
           read(66,*) jtestr(i)
        enddo
        close(66)

        OPEN(UNIT=66,FILE="../mathematica/jm1testi.txt",
     1      FORM="FORMATTED",STATUS="OLD",ACTION="READ")
        do i=1,ntest
           read(66,*) jtesti(i)
        enddo
        close(66)

        do i=1,ntest
           ztest(i)=rtest(i)*exp(ima*thetatest(i))
c           ztest(i)=rtest(i)
c           ztest(i) = ztest(i)/10
           write(79,*) dble(ztest(i)), imag(ztest(i))
        enddo

        do i=1,ntest
           jtest(i)=jtestr(i)+ima*jtesti(i)
        enddo

        ifexpon=0
        t1=second()
        do k = 1,10
        do i=1,ntest
           call abramm1(ztest(i),jtest2(i),ifexpon)
        enddo
        enddo
        t2=second()
        call prin2('evaluation time of J_{-1}=*',t2-t1,1)
        times(1,1)=t2-t1

        t1=second()
        do k = 1,10
        do i=1,ntest
           jtest3(i)=exp(ztest(i))
        enddo
        enddo
        t2=second()
        write(11,*) jtest3(1)
        call prin2('evaluation time of complex exponentials=*',t2-t1,1)
        times(1,2)=t2-t1

        dmax=0
        imax=0
        do i=1,ntest
           dd=abs((jtest(i)-jtest2(i))/jtest(i))
           if (dd .gt. 1d-14) then
              print *, dd, ztest(i)
           endif
           az=abs(ztest(i))
           if (az .lt. 1) then
              if (rerrors(1,1) .lt. dd) then
                 rerrors(1,1) = dd
              endif
           elseif (az .ge. 1 .and. az .lt. 3) then
              if (rerrors(2,1) .lt. dd) then
                 rerrors(2,1) = dd
              endif
           elseif (az .ge. 3 .and. az .lt. 15) then
              if (rerrors(3,1) .lt. dd) then
                 rerrors(3,1) = dd
              endif
           elseif (az .ge. 15 .and. az .lt. 120) then
              if (rerrors(4,1) .lt. dd) then
                 rerrors(4,1) = dd
              endif
           elseif (az .gt. 120) then
              if (rerrors(5,1) .lt. dd) then
                 rerrors(5,1) = dd
              endif
           endif

           if (dmax .lt. dd) then
              dmax=dd
              imax=i
           endif
        enddo

        print *, 'maximum relative error = ', dmax
        print *, 'occured at ', ztest(imax), imax
        print *, 'true value = ', jtest(imax)
        print *, 'computed value = ', jtest2(imax)
        print *, 'done with testing J_{-1}'

ccc        ntest=100 000
        OPEN(UNIT=66,FILE="../mathematica/r0test.txt",FORM="FORMATTED",
     1      STATUS="OLD",ACTION="READ")
        do i=1,ntest
           read(66,*) rtest(i)
        enddo
        close(66)

        OPEN(UNIT=66,FILE="../mathematica/theta0test.txt",
     1      FORM="FORMATTED",STATUS="OLD",ACTION="READ")
        do i=1,ntest
           read(66,*) thetatest(i)
        enddo
        close(66)

        OPEN(UNIT=66,FILE="../mathematica/j0testr.txt",FORM="FORMATTED",
     1      STATUS="OLD",ACTION="READ")
        do i=1,ntest
           read(66,*) jtestr(i)
        enddo
        close(66)

        OPEN(UNIT=66,FILE="../mathematica/j0testi.txt",FORM="FORMATTED",
     1      STATUS="OLD",ACTION="READ")
        do i=1,ntest
           read(66,*) jtesti(i)
        enddo
        close(66)

        do i=1,ntest
           ztest(i)=rtest(i)*exp(ima*thetatest(i))
c           ztest(i)=rtest(i)
c           ztest(i) = ztest(i)/10
           write(80,*) dble(ztest(i)), imag(ztest(i))
        enddo

        do i=1,ntest
           jtest(i)=jtestr(i)+ima*jtesti(i)
        enddo

        ifexpon=0
        t1=second()
        do k=1,10
        do i=1,ntest
           call abram0(ztest(i),jtest2(i),ifexpon)
        enddo
        enddo
        t2=second()
        call prin2('evaluation time of J_0=*',t2-t1,1)
        times(2,1)=t2-t1

        t1=second()
        do k=1,10
        do i=1,ntest
           jtest3(i)=exp(ztest(i))
        enddo
        enddo
        t2=second()
        write(11,*) jtest3(1)
        call prin2('evaluation time of complex exponentials=*',t2-t1,1)
        times(2,2)=t2-t1

        dmax=0
        imax=0
        do i=1,ntest
           dd=abs((jtest(i)-jtest2(i))/jtest(i))
           if (dd .gt. 1d-14) then
              print *, dd, ztest(i)
           endif
           az=abs(ztest(i))
           if (az .lt. 1) then
              if (rerrors(1,2) .lt. dd) then
                 rerrors(1,2) = dd
              endif
           elseif (az .ge. 1 .and. az .lt. 3) then
              if (rerrors(2,2) .lt. dd) then
                 rerrors(2,2) = dd
              endif
           elseif (az .ge. 3 .and. az .lt. 15) then
              if (rerrors(3,2) .lt. dd) then
                 rerrors(3,2) = dd
              endif
           elseif (az .ge. 15 .and. az .lt. 120) then
              if (rerrors(4,2) .lt. dd) then
                 rerrors(4,2) = dd
              endif
           elseif (az .gt. 120) then
              if (rerrors(5,2) .lt. dd) then
                 rerrors(5,2) = dd
              endif
           endif

           if (dmax .lt. dd) then
              dmax=dd
              imax=i
           endif
        enddo

        print *, 'maximum relative error = ', dmax
        print *, 'occured at ', ztest(imax), imax
        print *, 'true value = ', jtest(imax)
        print *, 'computed value = ', jtest2(imax)
        print *, 'done with testing J_0'

ccc        ntest=100 000
        OPEN(UNIT=66,FILE="../mathematica/r1test.txt",
     1      FORM="FORMATTED",STATUS="OLD",ACTION="READ")
        do i=1,ntest
           read(66,*) rtest(i)
        enddo
        close(66)

        OPEN(UNIT=66,FILE="../mathematica/theta1test.txt",
     1      FORM="FORMATTED",STATUS="OLD",ACTION="READ")
        do i=1,ntest
           read(66,*) thetatest(i)
        enddo
        close(66)

        OPEN(UNIT=66,FILE="../mathematica/j1testr.txt",
     1      FORM="FORMATTED",STATUS="OLD",ACTION="READ")
        do i=1,ntest
           read(66,*) jtestr(i)
        enddo
        close(66)

        OPEN(UNIT=66,FILE="../mathematica/j1testi.txt",
     1      FORM="FORMATTED",STATUS="OLD",ACTION="READ")
        do i=1,ntest
           read(66,*) jtesti(i)
        enddo
        close(66)

        do i=1,ntest
           ztest(i)=rtest(i)*exp(ima*thetatest(i))
c           ztest(i)=rtest(i)
c           ztest(i) = ztest(i)/10
           write(81,*) dble(ztest(i)), imag(ztest(i))
        enddo

        do i=1,ntest
           jtest(i)=jtestr(i)+ima*jtesti(i)
        enddo

        ifexpon=0
        t1=second()
        do k=1,10
        do i=1,ntest
           call abram1(ztest(i),jtest2(i),ifexpon)
        enddo
        enddo
        t2=second()
        call prin2('evaluation time of J_1=*',t2-t1,1)
        times(3,1)=t2-t1

        t1=second()
        do k=1,10
        do i=1,ntest
           jtest3(i)=exp(ztest(i))
        enddo
        enddo
        t2=second()
        write(11,*) jtest3(1)
        call prin2('evaluation time of complex exponentials=*',t2-t1,1)
        times(3,2)=t2-t1

        dmax=0
        imax=0
        do i=1,ntest
           dd=abs((jtest(i)-jtest2(i))/jtest(i))
           if (dd .gt. 1d-14) then
              print *, dd, ztest(i)
           endif
           az=abs(ztest(i))
           if (az .lt. 1) then
              if (rerrors(1,3) .lt. dd) then
                 rerrors(1,3) = dd
              endif
           elseif (az .ge. 1 .and. az .lt. 3) then
              if (rerrors(2,3) .lt. dd) then
                 rerrors(2,3) = dd
              endif
           elseif (az .ge. 3 .and. az .lt. 15) then
              if (rerrors(3,3) .lt. dd) then
                 rerrors(3,3) = dd
              endif
           elseif (az .ge. 15 .and. az .lt. 120) then
              if (rerrors(4,3) .lt. dd) then
                 rerrors(4,3) = dd
              endif
           elseif (az .gt. 120) then
              if (rerrors(5,3) .lt. dd) then
                 rerrors(5,3) = dd
              endif
           endif

           if (dmax .lt. dd) then
              dmax=dd
              imax=i
           endif
        enddo

        print *, 'maximum relative error = ', dmax
        print *, 'occured at ', ztest(imax), imax
        print *, 'true value = ', jtest(imax)
        print *, 'computed value = ', jtest2(imax)
        print *, 'done with testing J_1'

ccc        ntest=100 000
        OPEN(UNIT=66,FILE="../mathematica/r2test.txt",
     1      FORM="FORMATTED",STATUS="OLD",ACTION="READ")
        do i=1,ntest
           read(66,*) rtest(i)
        enddo
        close(66)

        OPEN(UNIT=66,FILE="../mathematica/theta2test.txt",
     1      FORM="FORMATTED",STATUS="OLD",ACTION="READ")
        do i=1,ntest
           read(66,*) thetatest(i)
        enddo
        close(66)

        OPEN(UNIT=66,FILE="../mathematica/j2testr.txt",
     1      FORM="FORMATTED",STATUS="OLD",ACTION="READ")
        do i=1,ntest
           read(66,*) jtestr(i)
        enddo
        close(66)

        OPEN(UNIT=66,FILE="../mathematica/j2testi.txt",
     1      FORM="FORMATTED",STATUS="OLD",ACTION="READ")
        do i=1,ntest
           read(66,*) jtesti(i)
        enddo
        close(66)

        do i=1,ntest
           ztest(i)=rtest(i)*exp(ima*thetatest(i))
c           ztest(i)=xtest(i)+ima*ytest(i)
c           ztest(i)=rtest(i)
c           ztest(i) = ztest(i)/10
           write(82,*) dble(ztest(i)), imag(ztest(i))
        enddo

        do i=1,ntest
           jtest(i)=jtestr(i)+ima*jtesti(i)
        enddo

        ifexpon=0
        t1=second()
        do k=1,10
        do i=1,ntest
           call abram2(ztest(i),jtest2(i),ifexpon)
        enddo
        enddo
        t2=second()
        call prin2('evaluation time of J_2=*',t2-t1,1)
        times(4,1)=t2-t1

        t1=second()
        do k=1,10
        do i=1,ntest
           jtest3(i)=exp(ztest(i))
        enddo
        enddo
        t2=second()
        write(11,*) jtest3(1)
        call prin2('evaluation time of complex exponentials=*',t2-t1,1)
        times(4,2)=t2-t1

        dmax=0
        imax=0
        do i=1,ntest
           dd=abs((jtest(i)-jtest2(i))/jtest(i))
           if (dd .gt. 1d-14) then
              print *, dd, ztest(i)
           endif
           az=abs(ztest(i))
           if (az .lt. 1) then
              if (rerrors(1,4) .lt. dd) then
                 rerrors(1,4) = dd
              endif
           elseif (az .ge. 1 .and. az .lt. 3) then
              if (rerrors(2,4) .lt. dd) then
                 rerrors(2,4) = dd
              endif
           elseif (az .ge. 3 .and. az .lt. 15) then
              if (rerrors(3,4) .lt. dd) then
                 rerrors(3,4) = dd
              endif
           elseif (az .ge. 15 .and. az .lt. 120) then
              if (rerrors(4,4) .lt. dd) then
                 rerrors(4,4) = dd
              endif
           elseif (az .gt. 120) then
              if (rerrors(5,4) .lt. dd) then
                 rerrors(5,4) = dd
              endif
           endif

           if (dmax .lt. dd) then
              dmax=dd
              imax=i
           endif
        enddo

        print *, 'maximum relative error = ', dmax
        print *, 'occured at ', ztest(imax), imax
        print *, 'true value = ', jtest(imax)
        print *, 'computed value = ', jtest2(imax)

        print *, 'done with testing J_2'

        write(29,*) '  &     S     &    Q_1 &     Q_2 &    Q_3 &    A\\'
 1600   format('J_{',I2,'}',2x,'&',4(2x,ES8.2,'&'), 2x, ES8.2, '\\')
        do i=1,4
           write(29,1600) i-2, (rerrors(j,i), j=1,5)
        enddo

        write(30,*) '  &     J_{-1} &     J_0 &     J_1 &     J_2\\'
 1800 format('T_E &', 3(2x,ES8.2,'&'),2x,ES8.2,'\\')
        do i=1,2
           write(30,1800) (times(j,i), j=1,4)
        enddo

        stop
        end

max_stars_repo_name: zgimbutas/abramowitz

max_issues_repo_name: zgimbutas/abramowitz

max_forks_repo_name: zgimbutas/abramowitz