#include "eiscor.h"
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! DOFDFI (Double Orthogonal hessenberg Factored Doubleshift Francis Iteration)
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! This routine computes one iteration of Francis' doubleshift
! algorithm on an orthogonal upper hessenberg matrix that is stored as a
! product of givens rotations and a diagonal matrix.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! INPUT VARIABLES:
!
!  COMPZ           CHARACTER
!                    'N': no eigenvectors
!                    'I': eigenvectors, initializing Z to the identity
!                    'V': eigenvectors, assume Z already initialized
!
!  N               INTEGER
!                    dimension of matrix
!
!  STR             INTEGER
!                    index of the top most givens rotation where
!                    the iteration begins
!
!  STP             INTEGER
!                    index of the bottom most givens rotation where
!                    the iteration ends
!
!  Q               REAL(8) array of dimension (2*(N-1))
!                    array of generators for givens rotations
!
!  D               REAL(8) array of dimension (2*N)
!                    array of generators for complex diagonal matrix
!
!  Z               REAL(8) array of dimension (N,N)
!                    if COMPZ = 'N' unused
!                    if COMPZ = 'I' stores eigenvectors in Z
!                    if COMPZ = 'V' update Z to store eigenvectors
!
!  ITCNT           INTEGER
!                   Contains the number of iterations since last deflation
!
! OUTPUT VARIABLES:
!
!  INFO            INTEGER
!                    INFO = 0 implies successful computation
!                    INFO = -1 implies COMPZ is invalid
!                    INFO = -3 implies STR is invalid
!                    INFO = -4 implies STP is invalid
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine DOFDFI(COMPZ,N,STR,STP,Q,D,Z,ITCNT,INFO)

  implicit none

  ! input variables
  character, intent(in) :: COMPZ
  integer, intent(in) :: N, STR, STP
  integer, intent(inout) :: ITCNT, INFO
  real(8), intent(inout) :: Q(2*(N-1)), D(2*N), Z(N,N)

  ! compute variables
  integer :: ii, ind
  real(8) :: s1, s2
  real(8) :: b1(2), b2(2), b3(2), temp(2)
  real(8) :: block(2,2)
  complex(8) :: eigs(2), h(2,2)

  ! initialize INFO
  INFO = 0

  ! check input in debug mode
  if (DEBUG) then

    ! check COMPZ
    if ((COMPZ.NE.'N').AND.(COMPZ.NE.'I').AND.(COMPZ.NE.'V')) then
      INFO = -1
      call UARERR(__FILE__,__LINE__,"COMPZ must be 'N', 'I' or 'V'",INFO,INFO)
      return
    end if

    ! check STR
    if ((STR < 1).OR.(STR > N-1)) then
      INFO = -3
      call UARERR(__FILE__,__LINE__,"STR must be 1 <= STR <= N-1",INFO,INFO)
      return
    end if

    ! check STP
    if ((STP < STR).OR.(STP > N-1)) then
      INFO = -4
      call UARERR(__FILE__,__LINE__,"STP must be STR <= STP <= N-1",INFO,INFO)
      return
    end if

  end if

  ! compute a nonzero shift
  ! shifts +1
  if((mod(ITCNT+1,11) == 0).AND.(ITCNT.NE.0))then
    eigs(1) = cmplx(1d0,0d0,kind=8)
    eigs(2) = cmplx(1d0,0d0,kind=8)

  ! shifts -1
  else if((mod(ITCNT+1,16) == 0).AND.(ITCNT.NE.0))then
    eigs(1) = cmplx(-1d0,0d0,kind=8)
    eigs(2) = cmplx(-1d0,0d0,kind=8)

  ! random shift
  else if((mod(ITCNT+1,21) == 0).AND.(ITCNT.NE.0))then
    call random_number(s1)
    call random_number(s2)
    eigs(1) = cmplx(s1,s2,kind=8)
    eigs(2) = cmplx(s1,-s2,kind=8)

  ! wilkinson shifts
  else
    ! get 2x2 block
    call DOFTDB(N,STP,Q,D,block,INFO)

    ! check INFO in debug mode
    if (DEBUG) then
      call UARERR(__FILE__,__LINE__,"DOFTDB failed",INFO,INFO)
      if (INFO.NE.0) then
        return
      end if
    end if

    ! compute eigenvalues and eigenvectors
    call DTTEEV(block,eigs,h,INFO)

    ! check INFO in debug mode
    if (DEBUG) then
      call UARERR(__FILE__,__LINE__,"DTTEEV failed",INFO,INFO)
      if (INFO.NE.0) then
        return
      end if
    end if

  end if

  ! two real shifts
  if ((aimag(eigs(1)).EQ.0d0).AND.(aimag(eigs(2)).EQ.0d0)) then

    ! first shift
    eigs(1) = cmplx(sign(1d0,dble(eigs(1))),0d0,kind=8)

    ! build first bulge
    temp(1) = dble(eigs(1))
    temp(2) = aimag(eigs(1))
    call DOFCFT('S',N,STR,Q,D,temp,b1,b2,INFO)

    ! check INFO in debug mode
    if (DEBUG) then
      call UARERR(__FILE__,__LINE__,"DOFCFT failed",INFO,INFO)
      if (INFO.NE.0) then
        return
      end if
    end if

    ! fusion to initialize first bulge
    ind = 2*(STR-1)
    b3(1) = b1(1)
    b3(2) = -b1(2)
    call DARFGR('R',b3,Q((ind+1):(ind+2)),INFO)

    ! check INFO in debug mode
    if (DEBUG) then
      call UARERR(__FILE__,__LINE__,"DARFGR failed",INFO,INFO)
      if (INFO.NE.0) then
        return
      end if
    end if

    ! first bulge update eigenvectors
    if (COMPZ .NE. 'N')then
      h(1,1) = b1(1)
      h(2,1) = b1(2)
      h(1,2) = -h(2,1)
      h(2,2) = h(1,1)
      Z(:,STR:(STR+1)) = matmul(Z(:,STR:(STR+1)),h)
    end if

    ! first bulge through D
    call DARGTD('R',D((ind+1):(ind+4)),b1,INFO)

    ! check INFO in debug mode
    if (DEBUG) then
      call UARERR(__FILE__,__LINE__,"DARGTD failed",INFO,INFO)
      if (INFO.NE.0) then
        return
      end if
    end if

    ! first bulge through Q
    call DARGTO(Q((ind+1):(ind+2)),Q((ind+3):(ind+4)),b1,INFO)

    ! check INFO in debug mode
    if (DEBUG) then
      call UARERR(__FILE__,__LINE__,"DARGTO failed",INFO,INFO)
      if (INFO.NE.0) then
        return
      end if
    end if

    ! second shift
    eigs(2) = cmplx(sign(1d0,dble(eigs(2))),0d0,kind=8)

    ! build bulge
    temp(1) = dble(eigs(2))
    temp(2) = aimag(eigs(2))
    call DOFCFT('S',N,STR,Q,D,temp,b2,b3,INFO)

    ! check INFO in debug mode
    if (DEBUG) then
      call UARERR(__FILE__,__LINE__,"DOFCFT failed",INFO,INFO)
      if (INFO.NE.0) then
        return
      end if
    end if

    ! fusion to initialize second bulge
    b3(1) = b2(1)
    b3(2) = -b2(2)
    call DARFGR('R',b3,Q((ind+1):(ind+2)),INFO)

    ! check INFO in debug mode
    if (DEBUG) then
      call UARERR(__FILE__,__LINE__,"DARFGR failed",INFO,INFO)
      if (INFO.NE.0) then
        return
      end if
    end if

    ! main chasing loop
    do ii=STR,(STP-2)

      ! set ind
      ind = 2*(ii-1)

      ! first bulge update eigenvectors
      if (COMPZ .NE. 'N')then
        h(1,1) = b1(1)
        h(2,1) = b1(2)
        h(1,2) = -h(2,1)
        h(2,2) = h(1,1)
        Z(:,(ii+1):(ii+2)) = matmul(Z(:,(ii+1):(ii+2)),h)
      end if

      ! first bulge through D
      call DARGTD('R',D((ind+3):(ind+6)),b1,INFO)

      ! check INFO in debug mode
      if (DEBUG) then
        call UARERR(__FILE__,__LINE__,"DARGTD failed",INFO,INFO)
        if (INFO.NE.0) then
          return
        end if
      end if

      ! first bulge through Q
      call DARGTO(Q((ind+3):(ind+4)),Q((ind+5):(ind+6)),b1,INFO)

      ! check INFO in debug mode
      if (DEBUG) then
        call UARERR(__FILE__,__LINE__,"DARGTO failed",INFO,INFO)
        if (INFO.NE.0) then
          return
        end if
      end if

      ! second bulge update eigenvectors
      if (COMPZ .NE. 'N')then
        h(1,1) = b2(1)
        h(2,1) = b2(2)
        h(1,2) = -h(2,1)
        h(2,2) = h(1,1)
        Z(:,ii:(ii+1)) = matmul(Z(:,ii:(ii+1)),h)
      end if

      ! second bulge through D
      call DARGTD('R',D((ind+1):(ind+4)),b2,INFO)

      ! check INFO in debug mode
      if (DEBUG) then
        call UARERR(__FILE__,__LINE__,"DARGTD failed",INFO,INFO)
        if (INFO.NE.0) then
          return
        end if
      end if

      ! second bulge through Q
      call DARGTO(Q((ind+1):(ind+2)),Q((ind+3):(ind+4)),b2,INFO)

      ! check INFO in debug mode
      if (DEBUG) then
        call UARERR(__FILE__,__LINE__,"DARGTO failed",INFO,INFO)
        if (INFO.NE.0) then
          return
        end if
      end if

    end do

    ! set ind
    ind = 2*(stp-2)

    ! first bulge update eigenvectors
    if (COMPZ .NE. 'N')then
      h(1,1) = b1(1)
      h(2,1) = b1(2)
      h(1,2) = -h(2,1)
      h(2,2) = h(1,1)
      Z(:,STP:(STP+1)) = matmul(Z(:,STP:(STP+1)),h)
    end if

    ! first bulge through D
    call DARGTD('R',D((ind+3):(ind+6)),b1,INFO)

    ! check INFO in debug mode
    if (DEBUG) then
      call UARERR(__FILE__,__LINE__,"DARGTD failed",INFO,INFO)
      if (INFO.NE.0) then
        return
      end if
    end if

    ! first bulge fuse with Q
    call DARFGR('L',Q((ind+3):(ind+4)),b1,INFO)

    ! check INFO in debug mode
    if (DEBUG) then
      call UARERR(__FILE__,__LINE__,"DARFGR failed",INFO,INFO)
      if (INFO.NE.0) then
        return
      end if
    end if

    ! second bulge update eigenvectors
    if (COMPZ .NE. 'N')then
      h(1,1) = b2(1)
      h(2,1) = b2(2)
      h(1,2) = -h(2,1)
      h(2,2) = h(1,1)
      Z(:,(STP-1):STP) = matmul(Z(:,(STP-1):STP),h)
    end if

    ! second bulge through D
    call DARGTD('R',D((ind+1):(ind+4)),b2,INFO)

    ! check INFO in debug mode
    if (DEBUG) then
      call UARERR(__FILE__,__LINE__,"DARGTD failed",INFO,INFO)
      if (INFO.NE.0) then
        return
      end if
    end if

    ! second bulge through Q
    call DARGTO(Q((ind+1):(ind+2)),Q((ind+3):(ind+4)),b2,INFO)

    ! check INFO in debug mode
    if (DEBUG) then
      call UARERR(__FILE__,__LINE__,"DARGTO failed",INFO,INFO)
      if (INFO.NE.0) then
        return
      end if
    end if

    ! set index
    ind = 2*(stp-1)

    ! last bulge update eigenvectors
    if (COMPZ .NE. 'N')then
      h(1,1) = b2(1)
      h(2,1) = b2(2)
      h(1,2) = -h(2,1)
      h(2,2) = h(1,1)
      Z(:,STP:(STP+1)) = matmul(Z(:,STP:(STP+1)),h)
    end if

    ! last bulge through D
    call DARGTD('R',D((ind+1):(ind+4)),b2,INFO)

    ! check INFO in debug mode
    if (DEBUG) then
      call UARERR(__FILE__,__LINE__,"DARGTD failed",INFO,INFO)
      if (INFO.NE.0) then
        return
      end if
    end if

    ! last bulge fuse with Q
    call DARFGR('L',Q((ind+1):(ind+2)),b2,INFO)

    ! check INFO in debug mode
    if (DEBUG) then
      call UARERR(__FILE__,__LINE__,"DARFGR failed",INFO,INFO)
      if (INFO.NE.0) then
        return
      end if
    end if

    ! update ITCNT
    ITCNT = ITCNT + 1

  ! complex conjugate pair
  else

    ! normalize shifts
    if (abs(eigs(1)).EQ.0d0) then
      eigs(1) = cmplx(0d0,1d0,kind=8)
    end if
    eigs(1) = eigs(1)/abs(eigs(1))

    ! build bulge
    temp(1) = dble(eigs(1))
    temp(2) = aimag(eigs(1))
    call DOFCFT('D',N,STR,Q,D,temp,b1,b2,INFO)

    ! check INFO in debug mode
    if (DEBUG) then
      call UARERR(__FILE__,__LINE__,"DOFCFT failed",INFO,INFO)
      if (INFO.NE.0) then
        return
      end if
    end if

    ! turnover to initialize bulge
    ind = 2*(STR-1)
    temp(1) = b2(1)
    temp(2) = -b2(2)
    b3(1) = b1(1)
    b3(2) = -b1(2)
    call DARGTO(temp,b3,Q((ind+1):(ind+2)),INFO)

    ! check INFO in debug mode
    if (DEBUG) then
      call UARERR(__FILE__,__LINE__,"DARGTO failed",INFO,INFO)
      if (INFO.NE.0) then
        return
      end if
    end if

    ! fusion to finish initialization
    call DARFGR('R',b3,Q((ind+3):(ind+4)),INFO)

    ! check INFO in debug mode
    if (DEBUG) then
      call UARERR(__FILE__,__LINE__,"DARFGR failed",INFO,INFO)
      if (INFO.NE.0) then
        return
      end if
    end if

    ! update b3
    b3 = Q((ind+1):(ind+2))

    ! update Q
    Q((ind+1):(ind+2)) = temp

    ! main chasing loop
    do ii=str,(stp-2)

      ! set ind
      ind = 2*(ii-1)

      ! first bulge update eigenvectors
      if (COMPZ .NE. 'N')then
        h(1,1) = b1(1)
        h(2,1) = b1(2)
        h(1,2) = -h(2,1)
        h(2,2) = h(1,1)
        Z(:,(ii+1):(ii+2)) = matmul(Z(:,(ii+1):(ii+2)),h)
      end if

      ! first bulge through D
      call DARGTD('R',D((ind+3):(ind+6)),b1,INFO)

      ! check INFO in debug mode
      if (DEBUG) then
        call UARERR(__FILE__,__LINE__,"DARGTD failed",INFO,INFO)
        if (INFO.NE.0) then
          return
        end if
      end if

      ! first bulge through Q
      call DARGTO(Q((ind+3):(ind+4)),Q((ind+5):(ind+6)),b1,INFO)

      ! check INFO in debug mode
      if (DEBUG) then
        call UARERR(__FILE__,__LINE__,"DARGTO failed",INFO,INFO)
        if (INFO.NE.0) then
          return
        end if
      end if

      ! second bulge update eigenvectors
      if (COMPZ .NE. 'N')then
        h(1,1) = b2(1)
        h(2,1) = b2(2)
        h(1,2) = -h(2,1)
        h(2,2) = h(1,1)
        Z(:,ii:(ii+1)) = matmul(Z(:,ii:(ii+1)),h)
      end if

      ! second bulge through D
      call DARGTD('R',D((ind+1):(ind+4)),b2,INFO)

      ! check INFO in debug mode
      if (DEBUG) then
        call UARERR(__FILE__,__LINE__,"DARGTD failed",INFO,INFO)
        if (INFO.NE.0) then
          return
        end if
      end if

      ! second bulge through Q
      call DARGTO(Q((ind+1):(ind+2)),Q((ind+3):(ind+4)),b2,INFO)

      ! check INFO in debug mode
      if (DEBUG) then
        call UARERR(__FILE__,__LINE__,"DARGTO failed",INFO,INFO)
        if (INFO.NE.0) then
          return
        end if
      end if

      ! push b3 down
      call DARGTO(b3,b1,b2,INFO)

      ! check INFO in debug mode
      if (DEBUG) then
        call UARERR(__FILE__,__LINE__,"DARGTO failed",INFO,INFO)
        if (INFO.NE.0) then
          return
        end if
      end if

    ! update bulges
      temp = b2
      b2 = b3
      b3 = b1
      b1 = temp

    end do

    ! set ind
    ind = 2*(stp-2)

    ! first bulge update eigenvectors
    if (COMPZ .NE. 'N')then
      h(1,1) = b1(1)
      h(2,1) = b1(2)
      h(1,2) = -h(2,1)
      h(2,2) = h(1,1)
      Z(:,STP:(STP+1)) = matmul(Z(:,STP:(STP+1)),h)
    end if

    ! first bulge through D
    call DARGTD('R',D((ind+3):(ind+6)),b1,INFO)

    ! check INFO in debug mode
    if (DEBUG) then
      call UARERR(__FILE__,__LINE__,"DARGTD failed",INFO,INFO)
      if (INFO.NE.0) then
        return
      end if
    end if

    ! first bulge fuse with Q
    call DARFGR('L',Q((ind+3):(ind+4)),b1,INFO)

    ! check INFO in debug mode
    if (DEBUG) then
      call UARERR(__FILE__,__LINE__,"DARFGR failed",INFO,INFO)
      if (INFO.NE.0) then
        return
      end if
    end if

    ! second bulge update eigenvectors
    if (COMPZ .NE. 'N')then
      h(1,1) = b2(1)
      h(2,1) = b2(2)
      h(1,2) = -h(2,1)
      h(2,2) = h(1,1)
      Z(:,(STP-1):STP) = matmul(Z(:,(STP-1):STP),h)
    end if

    ! second bulge through D
    call DARGTD('R',D((ind+1):(ind+4)),b2,INFO)

    ! check INFO in debug mode
    if (DEBUG) then
      call UARERR(__FILE__,__LINE__,"DARGTD failed",INFO,INFO)
      if (INFO.NE.0) then
        return
      end if
    end if

    ! second bulge through Q
    call DARGTO(Q((ind+1):(ind+2)),Q((ind+3):(ind+4)),b2,INFO)

    ! check INFO in debug mode
    if (DEBUG) then
      call UARERR(__FILE__,__LINE__,"DARGTO failed",INFO,INFO)
      if (INFO.NE.0) then
        return
      end if
    end if

    ! fuse b2 and b3
    call DARFGR('L',b3,b2,INFO)

    ! check INFO in debug mode
    if (DEBUG) then
      call UARERR(__FILE__,__LINE__,"DARFGR failed",INFO,INFO)
      if (INFO.NE.0) then
        return
      end if
    end if

    ! set index
    ind = 2*(stp-1)

    ! last bulge update eigenvectors
    if (COMPZ .NE. 'N')then
      h(1,1) = b3(1)
      h(2,1) = b3(2)
      h(1,2) = -h(2,1)
      h(2,2) = h(1,1)
      Z(:,STP:(STP+1)) = matmul(Z(:,STP:(STP+1)),h)
    end if

    ! last bulge through D
    call DARGTD('R',D((ind+1):(ind+4)),b3,INFO)

    ! check INFO in debug mode
    if (DEBUG) then
      call UARERR(__FILE__,__LINE__,"DARGTD failed",INFO,INFO)
      if (INFO.NE.0) then
        return
      end if
    end if

    ! last bulge fuse with Q
    call DARFGR('L',Q((ind+1):(ind+2)),b3,INFO)

    ! check INFO in debug mode
    if (DEBUG) then
      call UARERR(__FILE__,__LINE__,"DARFGR failed",INFO,INFO)
      if (INFO.NE.0) then
        return
      end if
    end if

    ! update ITCNT
    ITCNT = ITCNT + 1

  end if ! end of doubleshift step

end subroutine DOFDFI


text:
SUBROUTINE LUDCMPX (A , INDX , D , N , NP , ERROR )
IMPLICIT NONE
!PARAMETER DECLARATION
INTEGER, PARAMETER:: NMAX=50
!DATA DECLARATION
REAL(8) R0    /0.0D0/
REAL(8) R1    /1.0D0/
REAL(8) TINY  /1.0D-19/
CHARACTER	  NAME*6
DATA NAME/'LUDCMP'/
!SCALAR VARIABLES FROM ARGUMENTS
INTEGER N, NP
LOGICAL ERROR
REAL(8) D
!ARRAYS FROM ARGUMENTS
REAL(8), DIMENSION(NP,NP):: A
INTEGER, DIMENSION(N):: INDX
!LOCAL SCALAR VARIABLES
INTEGER I, J, K, IMAX
REAL(8) AAMAX, SUM, DUM
!LOCAL ARRAYS
REAL(8), DIMENSION(N):: VV
!INTIALIZE LOCAL SCALAR VARIABLES
I=0      ; J=0    ; K=0    ; IMAX=0
AAMAX=R0 ; SUM=R0 ; DUM=R0 ; VV=R0
!***********************************************************************
! Routine to do LU decomposition
! See NUMERICAL RECIPES p35.
!*ACRONYM
! LU_DeCOmPosition
!*DESCRIPTION
!*EXTERNAL
! Arrays
!=A      - LU decomposed matrix
!=INDX   - Permutation vector
! Variables
!=D      - Row interchange indicator
! N      - Size of the problem
! NP     - Physical size of A matrix ( NP >= N )
! ERROR  - Error flag
!*INTERNAL
! Arrays
! VV     - store the implicit scalling of each rows
!***********************************************************************
! Loop over each rows to get the implicit scalling factors
! --------------------------------------------------------
DO 12 I=1,N
	AAMAX=R0
        DO 11 J=1,N
        	IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
   	11   CONTINUE
        IF(AAMAX.EQ.R0)THEN
		ERROR=.TRUE.
          	GOTO 999
        ENDIF
        VV(I)=R1/AAMAX
12 CONTINUE
! Loop over columns of Crout's methods
! ------------------------------------
DO 19 J=1,N
	IF (J.GT.1) THEN
! Calculate upper triangle components Bij see equation(2.3.12), except i=j
        	DO 14 I=1,J-1
            		SUM=A(I,J)
            		IF (I.GT.1)THEN
              			DO 13 K=1,I-1
                			IF(ABS(A(I,K)).LT.TINY.AND.ABS(A(K,J)).LT.TINY)GOTO 13
! If A(I,K) and A(K,J) are very small, skip the following line
                			SUM=SUM-A(I,K)*A(K,J)
   				13         CONTINUE
              			A(I,J)=SUM
            		ENDIF
   		14     CONTINUE
	ENDIF
! Initialize for the search for largest pivot element
        AAMAX=R0
! Calculate diagonal components Bij see equation(2.3.12) ,i=j, and
! lower triangle components aij see equation 2.3.13, i=j+1,..N
        DO 16 I=J,N
        	SUM=A(I,J)
          	IF (J.GT.1)THEN
            		DO 15 K=1,J-1
              			IF(ABS(A(I,K)).LT.TINY.AND.ABS(A(K,J)).LT.TINY)GOTO 15
! If A(I,K) and A(K,J) are very small, skip the following line.
              			SUM=SUM-A(I,K)*A(K,J)
   			15       CONTINUE
            		A(I,J)=SUM
          	ENDIF
! Find largest pivot element
          	DUM=VV(I)*ABS(SUM)
          	IF (DUM.GE.AAMAX) THEN
            		IMAX=I
            		AAMAX=DUM
          	ENDIF
   	16   CONTINUE
! Interchange rows
        IF(J.NE.IMAX)THEN
        	DO 17 K=1,N
            		DUM=A(IMAX,K)
            		A(IMAX,K)=A(J,K)
            		A(J,K)=DUM
   		17     CONTINUE
! Change sign of row channge indicator and interchange scale factor
          	D=-D
          	VV(IMAX)=VV(J)
	ENDIF
! record interchange row number
        INDX(J)=IMAX
! Divided by the pivot element and get final aij values, see equation 2.3.13,
! i=j+1,..N
        IF(J.NE.N)THEN
        	IF(ABS(A(J,J)).LT.TINY)A(J,J)=TINY
          	DUM=R1/A(J,J)
          	DO 18 I=J+1,N
            		A(I,J)=A(I,J)*DUM
   		18     CONTINUE
	ENDIF
! Go back for next column
19 CONTINUE
IF(ABS(A(N,N)).LT.TINY)A(N,N)=TINY
999 CONTINUE
RETURN
END

max_stars_repo_name: iarlopes/GPROPT

max_issues_repo_name: iarlopes/GPROPT

max_forks_repo_name: iarlopes/GPROPT


text:
        PROGRAM MAIN

        IMPLICIT NONE
        REAL EPS
        EPS=1.0
        DO WHILE (1.0+EPS.GT.1.0)
            EPS=EPS/2.0
        END DO
        EPS=EPS*2

        PRINT *, EPS

        END PROGRAM MAIN

ext: f


text:
      MODULE m_kernel2
C.............................................................kernel2
c solution of 4 coupled equations (logariphmic mesh) for (l,\mu)
c 1->k1, 2->k2
c d/dxP1=[-k1]P1+exp(x)[e-v+2mc^2]/c^2Q1+exp(x)B/c^2\sigma_z[-1,-1]Q1
c d/dxQ1=[k1]Q1-exp(x)[e-v]P1+exp(x)B[\sigma_z[1,1]P1+\sigma_z[1,2]P2]
c d/dxP2=[-k2]P2+exp(x)[e-v+2mc^2]/c^2Q2+exp(x)B/c^2\sigma_z[-1,-1]Q2
c d/dxQ2=[k2]Q2-exp(x)[e-v]P2+exp(x)B[\sigma_z[2,1]P1+\sigma_z[2,2]P2]
c notations:
c Pi=exp(x)*gi
c Qi=exp(x)*cfi
c \gamma=sqrt(k^2-(Z/c)^2)
c v=0.5[V(1/2)+V(-1/2)]
c b=0.5[V(1/2)-V(-1/2)]
c Rydberg units: in charge
c Hartree units: com.
c NSOL= 2 -  4 equations
c ------------                                     a. shick KFA 1996
      CONTAINS
      SUBROUTINE kernel2(mrad,nsol,xmj,k1,k2,xx1,xx2,e,v,b,ri,dx,
     +                   nmatch,nstart,dp,dq,wp,wq)

      USE m_constants, ONLY : c_light
      USE m_diff
      IMPLICIT NONE
C     ..
C     .. Scalar Arguments ..
      INTEGER, INTENT (IN) :: mrad
      REAL dx,e,xmj
      INTEGER k1,k2,nmatch,nsol,nstart
C     ..
C     .. Array Arguments ..
      REAL b(mrad),dp(2,2,mrad),dq(2,2,mrad),ri(mrad),v(mrad),
     +     wp(2,2,mrad),wq(2,2,mrad),xx1(4),xx2(4)
C     ..
C     .. Local Scalars ..
      REAL abz,bh,bmn,bmnp1,bn,bnp1,c1,c2,c3,cc,csq,dp11,dp12,
     +     dp13,dp14,dp21,dp22,dp23,dp24,dq11,dq12,dq13,dq14,dq21,dq22,
     +     dq23,dq24,dxd8,expdxh,gk,gk1,gk1k,gkk1,gmk,gmk1,p1c,p2c,
     +     q1c,q2c,r,vh,vmn,vmnp1,vn,vnp1,xk,xk1,xm
      INTEGER i,iny,ir,j,jri,n,nrk
C     ..
C     .. Local Arrays ..
      REAL bm(mrad),p1(mrad),p1s(mrad),p2(mrad),p2s(mrad),q1(mrad),
     +     q1s(mrad),q2(mrad),q2s(mrad),vm(mrad)
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC exp,real,sqrt
C     ..
C     .. Data statements ..
      DATA abz/0.5/
C     ..
      cc = c_light(2.0)
c
      DO ir = 1,mrad
         p1(ir) = 0.0
         q1(ir) = 0.0
         p1s(ir) = 0.0
         q1s(ir) = 0.0
c
         p2(ir) = 0.0
         q2(ir) = 0.0
         p2s(ir) = 0.0
         q2s(ir) = 0.0
c
         DO i = 1,2
            DO j = 1,2
c
               wp(i,j,ir) = 0.0
               wq(i,j,ir) = 0.0
c
               dp(i,j,ir) = 0.0
               dq(i,j,ir) = 0.0
            END DO
         END DO
      END DO
C
      csq = cc*cc
      expdxh = exp(dx/2.0)
      dxd8 = dx/8.0
      jri = nmatch
      nrk = jri
c-derivatives of potential
      CALL diff(
     >          mrad,v,dx,jri,
     <          vm)
      CALL diff(
     >          mrad,b,dx,jri,
     <          bm)
c- begin to solve coupled Dirac equation in "magnetic" field
      xk = real(k1)
      xk1 = real(k2)
      xm = xmj
      gk = -xm/ (xk+abz)
      gk1 = -xm/ (xk1+abz)
      gkk1 = -sqrt(1.0-gk*gk)
      gk1k = gkk1
      gmk = xm/ (xk-abz)
      gmk1 = xm/ (xk1-abz)
c
c     NB: V=R**2*POTENTIAL
c     NB: B=R**2*FIELD
c
      DO 30 iny = 1,nsol
C                                                                       WAB02140
         r = ri(nstart)
         vn = v(nstart)
         vmn = vm(nstart)
         bn = b(nstart)
         bmn = bm(nstart)
         IF (iny.EQ.1) THEN
            p1(nstart) = xx1(1)
            q1(nstart) = xx1(2)
            p2(nstart) = xx1(3)
            q2(nstart) = xx1(4)
         ELSE
            p1(nstart) = xx2(1)
            q1(nstart) = xx2(2)
            p2(nstart) = xx2(3)
            q2(nstart) = xx2(4)
         END IF
c---->f.po solution of Dirac eq.
         c1 = vn/r**2 - e
         c2 = 1.0 - c1/csq
c      C2=2.0-C1/CSQ !H
         c3 = bn/r/r
         p1s(nstart) = -xk*p1(nstart) + r* (c2+c3/csq*gmk)*q1(nstart)
         q1s(nstart) = xk*q1(nstart) + r*
     +                 ((c1+c3*gk)*p1(nstart)+c3*gkk1*p2(nstart))
C
         p2s(nstart) = -xk1*p2(nstart) + r* (c2+c3/csq*gmk1)*q2(nstart)
         q2s(nstart) = xk1*q2(nstart) +
     +                 r* ((c1+c3*gk1)*p2(nstart)+c3*gk1k*p1(nstart))
C   ******************************************************************
C   *                                                                *
C   *     THE RUNGE-KUTTA METHOD IS USED IN THE <NKR> FIRST STEPS    *
C   *                                                                *
C   ******************************************************************
C
         n = nstart
c-> y1
   10    p1c = p1(n)
         q1c = q1(n)
         p2c = p2(n)
         q2c = q2(n)
c-> k1=f(x1,y1)
         dp11 = dx* (-xk*p1c+r* (c2+c3/csq*gmk)*q1c)
         dq11 = dx* (xk*q1c+r* ((c1+c3*gk)*p1c+c3*gkk1*p2c))
         dp21 = dx* (-xk1*p2c+r* (c2+c3/csq*gmk1)*q2c)
         dq21 = dx* (xk1*q2c+r* ((c1+c3*gk1)*p2c+c3*gk1k*p1c))
c-> y11=y1+0.5*k1
         p1c = p1c + 0.5*dp11
         q1c = q1c + 0.5*dq11
         p2c = p2c + 0.5*dp21
         q2c = q2c + 0.5*dq21
c-> x1'=x1+0.5h
         r = r*expdxh
c-> interpolation of V and B for x1'
         vnp1 = v(n+1)
         vmnp1 = vm(n+1)
         bnp1 = b(n+1)
         bmnp1 = bm(n+1)
         vh = (vn+vnp1)*0.5 + (vmn-vmnp1)*dxd8
         bh = (bn+bnp1)*0.5 + (bmn-bmnp1)*dxd8
c-> k2=f(x1',y11)
         c1 = vh/r/r - e
         c2 = 1.0 - c1/csq
c      C2=2.0-C1/CSQ !H
         c3 = bh/r/r
C                                                                       WAB02850
         dp12 = dx* (-xk*p1c+r* (c2+c3/csq*gmk)*q1c)
         dq12 = dx* (xk*q1c+r* ((c1+c3*gk)*p1c+c3*gkk1*p2c))
         dp22 = dx* (-xk1*p2c+r* (c2+c3/csq*gmk1)*q2c)
         dq22 = dx* (xk1*q2c+r* ((c1+c3*gk1)*p2c+c3*gk1k*p1c))
c-> y12=y1+0.5*k2 (y12=y11-0.5*k1+0.5*k2)
         p1c = p1c + 0.5* (dp12-dp11)
         q1c = q1c + 0.5* (dq12-dq11)
         p2c = p2c + 0.5* (dp22-dp21)
         q2c = q2c + 0.5* (dq22-dq21)
c-> k3=f(x1',y12)
         dp13 = dx* (-xk*p1c+r* (c2+c3/csq*gmk)*q1c)
         dq13 = dx* (xk*q1c+r* ((c1+c3*gk)*p1c+c3*gkk1*p2c))
         dp23 = dx* (-xk1*p2c+r* (c2+c3/csq*gmk1)*q2c)
         dq23 = dx* (xk1*q2c+r* ((c1+c3*gk1)*p2c+c3*gk1k*p1c))
c-> y13=y1+k3 (y13=y12-0.5*k2+k3)
         p1c = p1c + dp13 - 0.5*dp12
         q1c = q1c - 0.5*dq12 + dq13
         p2c = p2c + dp23 - 0.5*dp22
         q2c = q2c - 0.5*dq22 + dq23
c-> x2=x1+h
         n = n + 1
         r = ri(n)
c-> k4=f(x2,y13)
         c1 = vnp1/r/r - e
         c2 = 1.0 - c1/csq
c      C2=2.0-C1/CSQ !H
         c3 = bnp1/r/r
c
         dp14 = dx* (-xk*p1c+r* (c2+c3/csq*gmk)*q1c)
         dq14 = dx* (xk*q1c+r* ((c1+c3*gk)*p1c+c3*gkk1*p2c))
         dp24 = dx* (-xk1*p2c+r* (c2+c3/csq*gmk1)*q2c)
         dq24 = dx* (xk1*q2c+r* ((c1+c3*gk1)*p2c+c3*gk1k*p1c))
c-> y2=y1+1/6*(k1+2*k2+2*k3+k4)
         p1(n) = p1(n-1) + (dp11+2.0* (dp12+dp13)+dp14)/6.0
         q1(n) = q1(n-1) + (dq11+2.0* (dq12+dq13)+dq14)/6.0
         p2(n) = p2(n-1) + (dp21+2.0* (dp22+dp23)+dp24)/6.0
         q2(n) = q2(n-1) + (dq21+2.0* (dq22+dq23)+dq24)/6.0
c-> f(x2,y2)
         p1s(n) = -xk*p1(n) + r* (c2+c3/csq*gmk)*q1(n)
         q1s(n) = xk*q1(n) + r* ((c1+c3*gk)*p1(n)+c3*gkk1*p2(n))
         p2s(n) = -xk1*p2(n) + r* (c2+c3/csq*gmk1)*q2(n)
         q2s(n) = xk1*q2(n) + r* ((c1+c3*gk1)*p2(n)+c3*gk1k*p1(n))
c-> redefinition of V ind B (2->1)
         vn = vnp1
         bn = bnp1
         vmn = vmnp1
         bmn = bmnp1
C                                                                       WAB03340
         IF (n-nrk.LT.0) GOTO 10
c---->Milne's method off now!
c
c storage of wave functions and their derivatives.
         DO ir = 1,jri
            wp(1,iny,ir) = p1(ir)
            wp(2,iny,ir) = p2(ir)
            wq(1,iny,ir) = q1(ir)
            wq(2,iny,ir) = q2(ir)
c derivatives
            dp(1,iny,ir) = p1s(ir)
            dp(2,iny,ir) = p2s(ir)
            dq(1,iny,ir) = q1s(ir)
            dq(2,iny,ir) = q2s(ir)
         END DO
C***********************************************************            WAB05210
   30 CONTINUE
c      write(oUnit,*) wp

      END SUBROUTINE kernel2
      END MODULE m_kernel2

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text:
      SUBROUTINE setsd(iseed1,iseed2)
C**********************************************************************
C
C     SUBROUTINE SETSD(ISEED1,ISEED2)
C               SET S-ee-D of current generator
C
C     Resets the initial  seed of  the current  generator to  ISEED1 and
C     ISEED2. The seeds of the other generators remain unchanged.
C
C     This is a transcription from Pascal to Fortran of routine
C     Set_Seed from the paper
C
C     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
C     with Splitting Facilities." ACM Transactions on Mathematical
C     Software, 17:98-111 (1991)
C
C
C                              Arguments
C
C
C     ISEED1 -> First integer seed
C                                   INTEGER ISEED1
C
C     ISEED2 -> Second integer seed
C                                   INTEGER ISEED1
C
C**********************************************************************
C     .. Parameters ..
      INTEGER numg
      PARAMETER (numg=32)
C     ..
C     .. Scalar Arguments ..
      INTEGER iseed1,iseed2
C     ..
C     .. Scalars in Common ..
      INTEGER a1,a1vw,a1w,a2,a2vw,a2w,m1,m2
C     ..
C     .. Arrays in Common ..
      INTEGER cg1(numg),cg2(numg),ig1(numg),ig2(numg),lg1(numg),
     +        lg2(numg)
      LOGICAL qanti(numg)
C     ..
C     .. Local Scalars ..
      INTEGER g
C     ..
C     .. External Functions ..
      LOGICAL qrgnin
      EXTERNAL qrgnin
C     ..
C     .. External Subroutines ..
      EXTERNAL getcgn,initgn
C     ..
C     .. Common blocks ..
      COMMON /globe/m1,m2,a1,a2,a1w,a2w,a1vw,a2vw,ig1,ig2,lg1,lg2,cg1,
     +       cg2,qanti
C     ..
C     .. Save statement ..
      SAVE /globe/
C     ..
C     .. Executable Statements ..
C     Abort unless random number generator initialized
      IF (qrgnin()) GO TO 10
      WRITE (*,*) ' SETSD called before random number generator ',
     +  ' initialized -- abort!'
      STOP ' SETSD called before random number generator initialized'

   10 CALL getcgn(g)
      ig1(g) = iseed1
      ig2(g) = iseed2
      CALL initgn(-1)
      RETURN

      END

max_stars_repo_name: ppernot/MC-ChemDB

max_issues_repo_name: ppernot/MC-ChemDB

max_forks_repo_name: ppernot/MC-ChemDB


text:
      SUBROUTINE version
      COMMON /ver/ver
      CHARACTER*80 ver
      DATA ver
     &/' SIOSEIS ver 2020.5.1 (17 Nov. 2020)'/
      PRINT *,ver
      RETURN
      END


! chk_inquire_size.f90 --
!     Check: does the compiler support the SIZE= keyword in the INQUIRE statement?
!
program chk_size
    implicit none

    integer :: filesize
    logical :: exists

    write( *, '(a)' ) 'Examine the size of the source file:'

    inquire( file = 'chk_inquire_size.f90', size = filesize, exist = exists )

    if ( exists ) then
        write( *, '(a,i0)') '    Reported file size: ', filesize
        if ( filesize < 0 ) then
            write( *, '(a,i0)') '    Note: there is something wrong - the file exists and'
            write( *, '(a,i0)') '          should therefore be at least 0 bytes in size'
        endif
    else
        write( *, '(a)')   '    The file "chk_inquire_size,f90" does not seem to exist'
        write( *, '(a)')   '    Please check this'
    endif

    open( 10, file = 'chk_inquire_size.f90' )
    inquire( 10, size = filesize, exist = exists )

    write( *, '(a)' ) 'After opening the file:'
    if ( exists ) then
        write( *, '(a,i0)') '    Reported file size: ', filesize
        if ( filesize < 0 ) then
            write( *, '(a,i0)') '    Note: there is something wrong - the file exists and'
            write( *, '(a,i0)') '          should therefore be at least 0 bytes in size'
        endif
    else
        write( *, '(a)')   '    The file "chk_inquire_size,f90" is reported as non-existent'
        write( *, '(a)')   '    There is something seriously wrong'
    endif
end program chk_size


text:
	SUBROUTINE LC_COUN  ( stcn, cnflag, iret )
C************************************************************************
C* LC_COUN								*
C*									*
C* This subroutine checks STCN to see if it is a country abbreviation.	*
C* The following countries are currently recognized:			*
C*									*
C*	US	United States		CN	Canada			*
C*	MX	Mexico			BW	Bangladesh		*
C*	AU	Australia		CI	China			*
C*									*
C* Countries whose abbreviations will conflict with US state names	*
C* should not be added to this list.					*
C*									*
C* LC_COUN  ( STCN, CNFLAG, IRET )					*
C*									*
C* Input parameters:							*
C*	STCN		CHAR*		State / country abbreviation	*
C*									*
C* Output parameters:							*
C*	CNFLAG		LOGICAL		Country flag			*
C*	IRET		INTEGER		Return code			*
C*				    	   0 = normal return		*
C**									*
C* Log:									*
C* I. Graffman/RDS	 7/87						*
C* M. desJardins/GSFC	 6/88	Cleaned up				*
C************************************************************************
	PARAMETER	( NCOUN = 6 )
	CHARACTER*(*)	stcn
	LOGICAL		cnflag
C*
	CHARACTER	cnlist (NCOUN)*12
	DATA cnlist	/ 'US', 'CN',' MX', 'BW', 'CI', 'AU' /
C------------------------------------------------------------------------
	iret   = 0
	cnflag = .false.
	icn    = 0
C
C*	Check to see if this is a country.
C
	DO  i = 1, NCOUN
	    IF  ( stcn .eq. cnlist (i) )  cnflag = .true.
	END DO
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   FITNESS ROUTINES   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! These are functions to load the data from the GMM and then to calc         !!
!! the gradient and the energy at these positions.                            !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

module fitness
    integer, parameter :: NDIMENSION = 7, NCOMPONENT = 8
    integer :: ORIG_DIM = 7
    integer :: SCALING_METHOD = 2
    ! 0 is no scaling, 1 is the symmetric
    ! 2 is the general scaling
    ! ONLY NEED ORIG DIM IF SCALING METHOD IS 1
    ! TODO: This definitely needs neatening up at some point!

    double precision, allocatable :: WEIGHT(:), DETERMINANTS(:), MEANVECTOR(:,:)
    double precision, allocatable :: MATRIX(:,:,:), SIG(:), SCALING(:)
    real, parameter :: PI = 3.1415927

    CONTAINS

    subroutine fitness_init()

        INTEGER J1, J2, J3
        DOUBLE PRECISION PI
        character(len=50) :: filename
        INTEGER TEST(2)

        ALLOCATE(SIG(NCOMPONENT))
        ALLOCATE(WEIGHT(NCOMPONENT))
        ALLOCATE(DETERMINANTS(NCOMPONENT))
        ALLOCATE(MEANVECTOR(NDIMENSION,NCOMPONENT))
        ALLOCATE(MATRIX(NDIMENSION, NDIMENSION,NCOMPONENT))
        ALLOCATE(SCALING(NCOMPONENT))

        OPEN(UNIT=222,FILE='data/fortran/weight.out',status='old')
        READ(222,*) WEIGHT(1:NCOMPONENT)
        CLOSE(222)

        OPEN(UNIT=222,FILE='data/fortran/determinants.out',status='old')
        READ(222,*) DETERMINANTS(1:NCOMPONENT)
        CLOSE(222)

        OPEN(UNIT=222,FILE='data/fortran/meanvector.out',status='old')
        DO J1 = 1, NCOMPONENT
            READ(222,*) MEANVECTOR(1:NDIMENSION,J1)
        END DO
        CLOSE(222)

        DO J2=1, NCOMPONENT

            write (filename, fmt='(A I8 A)') "data/fortran/matrix", J2, ".out"
            call strip_spaces(filename)
            filename = trim(filename)
            OPEN(UNIT=222,FILE=filename,status='old')

            DO J1 = 1, NDIMENSION
                READ(222,*) MATRIX(1:NDIMENSION,J1,J2)
            END DO
            CLOSE(222)
        END DO

        ! ONLY LOAD SIGMA VALUES IF WE NEED THEM
        IF(SCALING_METHOD == 1) THEN
            OPEN(UNIT=222,FILE='data/fortran/sig.out',status='old')
            READ(222,*) SIG(1:NCOMPONENT)
            CLOSE(222)
        END IF

        ! ONLY LOAD THE SCALINGS IF WE NEED THEM
        IF(SCALING_METHOD == 2) THEN
            OPEN(UNIT=222,FILE='data/fortran/scalings.out',status='old')
            READ(222,*) SCALING(1:NCOMPONENT)
            CLOSE(222)
        END IF

    end subroutine fitness_init

    ! Calculate the gradient of the fitness surface
    subroutine fitness_gradient(COORDS, E, GRAD)

        double precision, intent(in) ::  COORDS(NDIMENSION)
        double precision, intent(out) :: GRAD(NDIMENSION)
        double precision, intent(out) :: E
        INTEGER J1, J2, J3
        DOUBLE PRECISION PROD(NCOMPONENT), PROB(NCOMPONENT), VECT(NDIMENSION)

        E = 0.0D0
        PROD(:) = 0.0D0
        PROB(:) = 0.0D0
        GRAD = 0d0

        DO J1 = 1, NCOMPONENT
             ! Starting vector = 0 for this component
         VECT(:) = 0.0D0

         DO J2 = 1, NDIMENSION
            DO J3 = 1, NDIMENSION
                         ! This is doing the sigma * (x - mu)
               VECT(J2) = VECT(J2) + MATRIX(J2,J3, J1) * ( COORDS(J3)-MEANVECTOR(J3, J1) )
            ENDDO
                    ! This is doing th (x - mu) * that
            PROD(J1) = PROD(J1) + VECT(J2) * (COORDS(J2)-MEANVECTOR(J2,J1))
         ENDDO
             ! This calculates the probability and the energy

         IF (SCALING_METHOD == 0) THEN
                PROB(J1) = EXP( -PROD(J1)/2.0 ) / SQRT( DETERMINANTS(J1) )
             ELSE IF (SCALING_METHOD == 1) THEN
                PROB(J1) = EXP( -PROD(J1)/2.0 ) / SQRT( DETERMINANTS(J1) )*(1/( SQRT( 2*PI*SIG(J1) ) ))**(ORIG_DIM-NDIMENSION)
             ELSE
                  PROB(J1) = EXP( -PROD(J1)/2.0 ) / SQRT( DETERMINANTS(J1) )*SCALING(J1)
         END IF

         E = E + WEIGHT(J1)*PROB(J1)

             ! Once the probability has been calculated we can calculate the gradient for each direction
         DO J2 = 1, NDIMENSION
            GRAD(J2) = GRAD(J2) - PROB(J1) * VECT(J2) * WEIGHT(J1)
         END DO

        END  DO

        ! Logging the energy
        GRAD = -GRAD/E
        E = -log(E)

    end subroutine fitness_gradient

    ! This is a dirty subroutine for removing spaces from a string
    ! this is neccesary for the file names matrix1.out
    ! otherwise we end up with matrix   1.out
    subroutine strip_spaces(string)
        character(len=*) :: string
        integer :: string_len
        integer :: last, actual

        string_len = len (string)
        last = 1
        actual = 1

        do while (actual < string_len)
            if (string(last:last) == ' ') then
                actual = actual + 1
                string(last:last) = string(actual:actual)
                string(actual:actual) = ' '
            else
                last = last + 1
                if (actual < last) &
                    actual = last
            endif
        end do
    end subroutine strip_spaces

end module fitness

max_stars_repo_name: philip-pearce/learning-dynamical

max_issues_repo_name: philip-pearce/learning-dynamical

max_forks_repo_name: philip-pearce/learning-dynamical


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

20090605 18:46:58 nbsp Welcome to the Wiki. Do you need to recover the password to your Users/lreyes account? Users/JasonAller

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
! Rice test
! If a pointer component of a derived type makes a forward reference
! to another derived type, it is unparsed with type 'integer'.

program forward_ref_derived_type

   type :: t1
      type(t2), pointer :: c  ! type of 'c' becomes 'integer'
      type(t2), pointer :: d  !!!!
      type(t99), pointer :: x99
      type(t1), pointer :: y77
   end type

   type :: t2
      real ::  x
      type(t99), pointer :: z99
      type(t1) :: wxyz
   end type

   type :: t99
   	  type(t2) :: a
   	  type(t1) :: b
   	  type(t99), pointer :: p99
   end type

   type(t1) a
   a%c%x = 1.0

end program

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
C    @(#)swpjdx.f	20.3 2/13/96
        subroutine swpjdx(m,n)
C
      include 'ipfinc/parametr.inc'
      include 'ipfinc/ikk.inc'
C
        do 100 i=1,3
           itemp=jndx(i,m)
           jndx(i,m)=jndx(i,n)
           jndx(i,n)=itemp
  100   continue
        return
        end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
      SUBROUTINE NON_LINEAR13 (IS, M, U, X, Y, NQ)

      IMPLICIT REAL*8 (A-H, O-Z)
      PARAMETER (PI = 3.141592653589793238D0)
      COMMON /COUNT/ KTTRANS(256)
      DIMENSION U(1), X(1), Y(1)

      if(x.gt.0) then
         is = -1
      else
         is = 1
      endif

      IF (IS .NE. 1 .AND. IS .NE. -1)
     $  THEN
         print *, is
      ENDIF
      END


text:
C$Procedure      ZZEKTR32 ( EK tree, 3-2 merge )

      SUBROUTINE ZZEKTR32 ( HANDLE,  TREE,    LEFT,    MIDDLE,
     .                      RIGHT,   PARENT,  LPKIDX,  UNDRFL )

C$ Abstract
C
C     Execute a 3-2 merge:  merge three neighboring sibling nodes, two
C     of which contain the minimum number of keys and one of which
C     has an underflow of one key, into two nodes, each one
C     approximately full.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     EK
C
C$ Keywords
C
C     EK
C     PRIVATE
C
C$ Declarations

      INCLUDE 'ektree.inc'
      INCLUDE 'ekpage.inc'
      INCLUDE 'ektype.inc'

      INTEGER               HANDLE
      INTEGER               TREE
      INTEGER               LEFT
      INTEGER               MIDDLE
      INTEGER               RIGHT
      INTEGER               PARENT
      INTEGER               LPKIDX
      LOGICAL               UNDRFL

C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C     HANDLE     I   File handle.
C     TREE       I   Root of tree.
C     LEFT       I   Left sibling node.
C     MIDDLE     I   Middle sibling node.
C     RIGHT      I   Right sibling node.
C     PARENT     I   Common parent node.
C     LPKIDX     I   Node-relative index of left parent key of MIDDLE.
C     UNDRFL     O   Underflow flag.
C
C$ Detailed_Input
C
C     HANDLE         is a file handle of an EK open for write access.
C
C     TREE           is the root node number of the tree of interest.
C
C     LEFT           is the node number of the left node of a trio of
C                    siblings.  LEFT either contains the minimum
C                    allowed number of keys  or is underflowing by
C                    one key.
C
C     MIDDLE         is the node number of the middle node of a trio of
C                    siblings.  MIDDLE either contains the minimum
C                    allowed number of keys  or is underflowing by
C                    one key.
C
C     RIGHT          is the node number of the right node of a trio of
C                    siblings.  The total number of keys in nodes
C                    LEFT, MIDDLE and RIGHT amounts to an underflow of 1
C                    key.
C
C     PARENT         is the node number of the common parent of LEFT
C                    and RIGHT.
C
C     LPKIDX         is the node-relative index within PARENT of the
C                    left parent key of MIDDLE.  This key is the
C                    immediate predecessor of the first key in the
C                    subtree headed by MIDDLE.
C
C$ Detailed_Output
C
C     UNDRFL         is a logical flag indicating whether the parent
C                    node underflowed as a result of the 3-2 merge.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     1)  If HANDLE is invalid, the error will be diagnosed by routines
C         called by this routine.  The file will not be modified.
C
C     2)  If an I/O error occurs while reading the indicated file, the
C         error will be diagnosed by routines called by this routine.
C
C     3)  If LEFT and RIGHT are not neighboring siblings, the error
C         SPICE(BUG) is signalled.
C
C     4)  If either LEFT or RIGHT are not children of PARENT, the error
C         SPICE(BUG) is signalled.
C
C     5)  If the sum of the number of keys in LEFT and RIGHT does not
C         correspond to an underflow of exactly 1 key, the error
C         SPICE(BUG) is signalled.
C
C$ Files
C
C     See the EK Required Reading for a discussion of the EK file
C     format.
C
C$ Particulars
C
C     Deletions from an EK tree start at a leaf node.  If the node
C     underflows, the EK system attempts to shuffle keys at the leaf
C     level to resolve the underflow.  That attempt failing, the system
C     delegates the problem upward to the next higher level.
C
C     There are two ways to resolve underflow in a non-root node:
C     balance the underflowing node with one of its closest siblings
C     (see ZZEKTRBN), or if the closest siblings contain the minimum
C     number of keys, execute a 3-2 merge.
C
C     A 3-2 merge involves deletion of the middle node of a trio of
C     neighboring siblings, two of which contain the minimum
C     number of keys and one of which has an underflow of one key,
C     and redistributing the keys between the two remaining nodes
C     and their common parent so that each of the two remaining siblings
C     contains the maximum number of keys.  The parent loses a key in
C     the process.
C
C     After the 3-2 merge, the tree is balanced and the siblings
C     satisfy the key count invariants.  However, the parent of the
C     siblings may underflow by one key.
C
C     Below are the gory details concerning the actions of this routine.
C     All of the parameters referred to here (in capital letters) are
C     defined in the include file ektree.inc.
C
C     In a 3-2 merge:
C
C        - The left parent key of the middle child is rotated down
C          into the left child and is appended on the right to the key
C          set of that child.  The left parent key's data pointer moves
C          along with the key.  The child pointers of this parent key
C          do not move along with the key; these pointers point to the
C          left and middle child nodes.
C
C        - The keys of the middle child are divided into three sets:
C          a set to be rotated left through the parent node into the
C          left child, a singleton set consisting of a key to be moved
C          up into the parent, and a set of keys to be rotated right
C          through the parent into the right child.  The sizes of the
C          leftmost and rightmost of these sets differ by at most 1.
C
C        - The number of keys that are rotated left is picked so that
C          after the rotation, the size of the left node will be
C          (3*MNKEYC)/2.  The data pointers of these keys move along
C          with the keys.  All of the left and right child pointers
C          of these keys move along with the keys into the left child.
C
C        - The singleton key in the child moves up into the parent
C          node.  Its data pointer moves with it.  After the move into
C          the parent node, the left child pointer of this key points
C          to the left child; the right child pointer points to the
C          right child.
C
C        - The right parent key of the middle child is rotated right
C          into the right child.  The data pointer of this key moves
C          with the key.  The child pointers of this parent key
C          do not move along with the key; these pointers point to the
C          middle and right child nodes.
C
C        - The remaining keys in the middle child are rotated right
C          into the right child; these become the leftmost keys of
C          that child.  The data pointers of these keys move along
C          with them.  The child pointers of these keys also move
C          along with them.
C
C        - The right child ends up with
C
C             (3*MNKEYC) - (3*MNKEYC)/2
C
C          keys.  This may be deduced from the facts that the original
C          three children had between them an underflow of one key, the
C          parent lost a key, and the left child has (3*MNKEYC)/2 keys.
C
C          Since
C
C             MNKEYC    =  MNKIDC - 1
C                       =  (  ( 2*MXKIDC + 1 ) / 3   )  -  1
C                       =  ( 2*MXKIDC - 2 ) / 3
C                       =  ( 2*MXKEYC ) / 3
C
C          we have
C
C             3*MNKEYC  <  2*MXKEYC
C                       -
C
C          If 3*MNKEYC is odd, we have strict inequality and also
C
C
C             3*MNKEYC  =  2 * ( (3*MNKEYC)/2 )  +  1
C
C          so
C
C             3*MNKEYC  + 1  =  2 * ( (3*MNKEYC)/2 )  +  2
C
C                            =  2 * (  (3*MNKEYC)/2  +  1 )
C
C                            <  2 * MXKEYC
C                            -
C
C          So in this case, the larger of the child nodes, which has
C          size
C
C             (3*MNKEYC)/2  +  1
C
C          has a key count no greater than MXKEYC.
C
C          If 3*MNKEYC is even, then the left and right child are the
C          same size, and the inequality
C
C             3*MNKEYC  <  2*MXKEYC
C                       -
C
C          implies directly that both nodes have size no greater than
C          MXKEYC.
C
C          Since both child nodes have size approximately 3/2 * MNKEYC,
C          and since MNKEYC is approximately 40, the minimum size
C          constraints on the child nodes are easily met.
C
C
C     As the above description shows, the parent loses a key as a
C     result of a 3-2 merge.  This may cause the parent to underflow;
C     if it does, the underflow at the parent's level must be resolved.
C
C$ Examples
C
C     See ZZEKTRDL.
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     N.J. Bachman   (JPL)
C
C$ Version
C
C-    Beta Version 1.0.0, 16-NOV-1995 (NJB)
C
C-&

C
C     Non-SPICELIB functions
C
      INTEGER               ZZEKTRBS

C
C     Local parameters
C
      INTEGER               LSIZE
      PARAMETER           ( LSIZE  = (3*MNKEYC)/2 )

      INTEGER               RSIZE
      PARAMETER           ( RSIZE  =  3*MNKEYC - LSIZE )

C
C     Local variables
C
      INTEGER               C1PAGE ( PGSIZI )
      INTEGER               C2PAGE ( PGSIZI )
      INTEGER               C3PAGE ( PGSIZI )
      INTEGER               DATBAS
      INTEGER               I
      INTEGER               KEYBAS
      INTEGER               KIDBAS
      INTEGER               LEFTSZ
      INTEGER               LMIDSZ
      INTEGER               LPKEY
      INTEGER               LSIB
      INTEGER               MIDSIZ
      INTEGER               MSIB
      INTEGER               N
      INTEGER               NLKEYS
      INTEGER               NMKEYS
      INTEGER               NNODE
      INTEGER               NPKEYS
      INTEGER               NRKEYS
      INTEGER               PPAGE ( PGSIZI )
      INTEGER               PSIZE
      INTEGER               RBASE
      INTEGER               RMIDSZ
      INTEGER               ROOT
      INTEGER               RPKEY
      INTEGER               RSHIFT
      INTEGER               RSIB
      INTEGER               SIZBAS
      INTEGER               SUM

C
C     Use discovery check-in for speed.
C
C     The plan is to take three sibling nodes, two of which contain
C     the minimum number of keys and one of which is underflowing by one
C     key, and to merge these nodes into two nodes.  This process
C     reduces the number of nodes in the parent by one and may cause the
C     parent to underflow.
C
C     After the merge, the sum of the numbers of keys in the two
C     children will be exactly (3*MNKEYC).  The numbers of keys in the
C     left and right nodes will be, respectively:
C
C
C        LSIZE  =   INT (  (3*MNKEYC)/2  )
C        RSIZE  =   (3*MNKEYC)  - LSIZE
C
C     We need to be sure that LSIZE and RSIZE are in the range
C
C        MNKEYC : MXKEYC
C
C
C     The definition of LSIZE implies that
C
C        LSIZE  =  MNKEYC + INT ( MNKEYC/2 )
C
C
C     so
C
C        MNKEYC + INT ( MNKEYC/2 )  <  LSIZE <  (3/2)*MNKEYC
C                                   -        -
C
C     and since
C
C        MNKEYC  =  MNKIDC - 1
C                =  INT ( ( 2*MXKIDC + 1 ) / 3 ) - 1
C                =  INT ( ( 2*MXKEYC + 3 ) / 3 ) - 1
C                =  INT ( ( 2*MXKEYC     ) / 3 )
C
C     we have
C
C        (3/2) * MNKEYC  =  (3/2) * INT ( (2*MXKEYC) / 3 )  <  MXKEYC
C                                                           -
C
C     Thus  LSIZE is guaranteed to be in range.
C
C     When MNKEYC is even, RSIZE is equal to LSIZE and thus is
C     within bounds.  When MNKEYC is odd, RSIZE exceeds LSIZE by 1, so
C
C        MNKEYC  <  RSIZE
C
C
C     It remains to be shown that
C
C        RSIZE   <  MXKEYC
C                -
C
C     when MNKEYC is odd.  When this is the case, the quantity
C
C        (3/2) * MNKEYC
C
C     is not an integer and therefore is strictly less than MXKEYC.
C     This quantity is also greater than LSIZE, so we conclude that
C
C        LSIZE  <  MXKEYC - 1
C               -
C
C     Since RSIZE exceeds LSIZE by 1, we have
C
C        RSIZE  <  MXKEYC
C               -
C
C     as we claimed.
C
C
C     All right, read in the child and parent pages.
C
      CALL ZZEKPGRI ( HANDLE, LEFT,   C1PAGE )
      CALL ZZEKPGRI ( HANDLE, MIDDLE, C2PAGE )
      CALL ZZEKPGRI ( HANDLE, RIGHT,  C3PAGE )
      CALL ZZEKPGRI ( HANDLE, PARENT, PPAGE  )

C
C     The actual addresses in the parent node depend on whether the
C     parent is the root.  Compute the necessary bases to avoid a lot
C     of cases.
C
      ROOT  =  TREE

      IF ( PARENT .EQ. ROOT ) THEN
         KEYBAS  =  TRKEYR
         DATBAS  =  TRDATR
         KIDBAS  =  TRKIDR
         SIZBAS  =  TRNKR
      ELSE
         KEYBAS  =  TRKEYC
         DATBAS  =  TRDATC
         KIDBAS  =  TRKIDC
         SIZBAS  =  TRNKC
      END IF

C
C     Check the left parent key of the middle child.
C
      PSIZE =  PPAGE ( SIZBAS )

      IF (  ( LPKIDX .LT. 1 ) .OR. ( LPKIDX .GT. PSIZE-1 )  )  THEN

         CALL CHKIN  ( 'ZZEKTR32'                                )
         CALL SETMSG ( 'Left parent key of MIDDLE is out of '   //
     .                 'range.  Value is #; valid range is 1:#'  )
         CALL ERRINT ( '#',   LPKIDX                             )
         CALL ERRINT ( '#',   PSIZE-1                            )
         CALL SIGERR ( 'SPICE(BUG)'                              )
         CALL CHKOUT ( 'ZZEKTR32'                                )
         RETURN

      END IF

C
C     Retain the left and right parent key values of the middle child.
C
      LPKEY  =  PPAGE( KEYBAS + LPKIDX    )
      RPKEY  =  PPAGE( KEYBAS + LPKIDX + 1)

C
C     Verify that LEFT, MIDDLE, and RIGHT are siblings, and that PARENT
C     is their common parent.
C
      LSIB  =  PPAGE ( KIDBAS+LPKIDX   )
      MSIB  =  PPAGE ( KIDBAS+LPKIDX+1 )
      RSIB  =  PPAGE ( KIDBAS+LPKIDX+2 )

      IF (      ( LSIB .NE. LEFT   )
     .     .OR. ( MSIB .NE. MIDDLE )
     .     .OR. ( RSIB .NE. RIGHT  )  ) THEN

         CALL CHKIN  ( 'ZZEKTR32'                                      )
         CALL SETMSG ( 'LEFT, RIGHT, MIDDLE, PARENT, and PKIDX are '  //
     .                 'inconsistent. LEFT = #; MIDDLE = #; '         //
     .                 'RIGHT = #; PARENT = #; LPKIDX = #; LSIB '     //
     .                 'derived from PARENT = #; MSIB = #; RSIB = #.'  )
         CALL ERRINT ( '#',   LEFT                                     )
         CALL ERRINT ( '#',   MIDDLE                                   )
         CALL ERRINT ( '#',   RIGHT                                    )
         CALL ERRINT ( '#',   PARENT                                   )
         CALL ERRINT ( '#',   LPKIDX                                   )
         CALL ERRINT ( '#',   LSIB                                     )
         CALL ERRINT ( '#',   MSIB                                     )
         CALL ERRINT ( '#',   RSIB                                     )
         CALL SIGERR ( 'SPICE(BUG)'                                    )
         CALL CHKOUT ( 'ZZEKTR32'                                      )
         RETURN

      END IF

C
C     Get the number of keys in the parent.
C
      IF ( PARENT .EQ. ROOT ) THEN
         NPKEYS = PPAGE ( TRNKR )
      ELSE
         NPKEYS = PPAGE ( TRNKC )
      END IF

C
C     Get the number of keys in each child.
C
      NLKEYS  =  C1PAGE ( TRNKC )
      NMKEYS  =  C2PAGE ( TRNKC )
      NRKEYS  =  C3PAGE ( TRNKC )

      SUM     =  NLKEYS + NMKEYS + NRKEYS

C
C     The sum of the number of keys in the three input nodes must
C     sum exactly to value representing an underflow level of 1 key.
C
      IF (  SUM  .NE.  ( 3*MNKEYC - 1 )  ) THEN

         CALL CHKIN  ( 'ZZEKTR32'                                      )
         CALL SETMSG ( 'Number of keys in nodes LEFT = #; in MIDDLE ' //
     .                 '= #; in RIGHT = #; counts summing to # were ' //
     .                 'expected.'                                     )
         CALL ERRINT ( '#',  NLKEYS                                    )
         CALL ERRINT ( '#',  NMKEYS                                    )
         CALL ERRINT ( '#',  NRKEYS                                    )
         CALL ERRINT ( '#',  3*MNKEYC-1                                )
         CALL SIGERR ( 'SPICE(BUG)'                                    )
         CALL CHKOUT ( 'ZZEKTR32'                                      )
         RETURN

      END IF

C
C     We're set to carry out the merge.  Here's an overview of what
C     gets moved where.
C
C        The left parent key of the middle node moves into the left
C        node, at the end of the node.
C
C        The first N-1 keys and N child pointers of the middle node get
C        moved into the left node, where
C
C           N  =  LSIZE - ( 1 + NLKEYS ) + 1
C
C        The Nth key of the middle node moves into the parent,
C        replacing the left parent key of the middle node.
C
C        The right parent key of the middle node moves into the right
C        node, at the beginning of the node.
C
C        The keys from position N+1 onward in the middle node, as
C        well as all of the remaining child pointers, move into the
C        right node, at the beginning.
C
C        The right parent key's location is filled in by shifting
C        the keys, data pointers, and child pointers in the parent
C        to the left by one position.  The child pointer removed by this
C        operation is the pointer to the middle child.
C
C        The middle child node disappears.
C
C     Before re-arranging things, we'll need to have on hand the key
C     counts for various sets of keys.  We'll use the variable LEFTSZ
C     for the number of keys in the subtree headed by LEFT.  We'll
C     use the variable LMIDSZ to refer to the `subtree' headed by
C     the set of keys in the middle node that will be shifted into
C     the left child.  The variable RMSIZE will represent the size of
C     the key set moved from the middle child into the right child.
C     MIDSIZ will be the key count for the subtree headed by the middle
C     child.
C
C     Consistent with usage above, the variable N will represent
C     the index of the key in the middle node that will rapturously
C     ascend into the parent.
C

      IF ( LPKIDX .EQ. 1 ) THEN
         LEFTSZ  =  LPKEY - 1
      ELSE
         LEFTSZ  =  LPKEY - PPAGE(KEYBAS+LPKIDX-1) - 1
      END IF

      N       =  LSIZE  -  ( 1 + NLKEYS ) + 1

      LMIDSZ  =  C2PAGE(TRKEYC+N)  - 1
      MIDSIZ  =  RPKEY  -  LPKEY   - 1
      RMIDSZ  =  MIDSIZ -  LMIDSZ  - 1

C
C     Move the left parent key into the left child.  The key itself
C     doesn't really move; its value is simply re-assigned.  The
C     data pointer is copied, however.  The child pointer at location
C     LSIZE+1 is unaffected by this move.
C
      C1PAGE ( TRKEYC + NLKEYS + 1 )  =  LEFTSZ + 1
      C1PAGE ( TRDATC + NLKEYS + 1 )  =  PPAGE ( DATBAS + LPKIDX )

C
C     Move the first N-1 keys and data pointers, and the first N
C     child pointers, from the middle child into the left
C     child.  The moved keys will gain LEFTSZ + 1 predecessors.
C
      DO I = 1, N-1
         C1PAGE(TRKEYC+NLKEYS+1+I) =  C2PAGE(TRKEYC+I) +  LEFTSZ + 1
      END DO

      CALL MOVEI ( C2PAGE(TRDATC+1),  N-1,  C1PAGE(TRDATC+NLKEYS+2) )
      CALL MOVEI ( C2PAGE(TRKIDC+1),  N,    C1PAGE(TRKIDC+NLKEYS+2) )

C
C     Set the key count in the left child.
C
      C1PAGE ( TRNKC )  =  LSIZE

C
C     The left child is complete.  Now it's time to set up the right
C     child.  First off, we'll shift the node's contents to the right
C     by the number of new keys we're going to insert.  Shift the
C     rightmost elements first.  The shifted keys will gain RMIDSZ+1
C     predecessors, so  we adjust the keys as we shift them.
C
      RSHIFT  =  NMKEYS - N + 1

      DO I = NRKEYS, 1, -1
         C3PAGE(TRKEYC+RSHIFT+I) =  C3PAGE(TRKEYC+I) + RMIDSZ + 1
      END DO

      DO I = NRKEYS, 1, -1
         C3PAGE(TRDATC+RSHIFT+I) =  C3PAGE(TRDATC+I)
      END DO

      DO I = NRKEYS+1, 1, -1
         C3PAGE(TRKIDC+RSHIFT+I) =  C3PAGE(TRKIDC+I)
      END DO

C
C     The key at location RSHIFT receives the former right parent key
C     of the middle child.   The key value is simply assigned; the
C     data pointer is copied.  The child pointer at location RSHIFT
C     will be set later.
C
      C3PAGE ( TRKEYC + RSHIFT )  =  RMIDSZ + 1
      C3PAGE ( TRDATC + RSHIFT )  =  PPAGE ( DATBAS+LPKIDX+1 )

C
C     The first RSHIFT-1 locations in the right child are filled in
C     with data from the middle child.  The moved keys lose LMIDSZ+1
C     precedessors.
C
      DO I = 1, RSHIFT-1
         C3PAGE(TRKEYC+I) = C2PAGE(TRKEYC+N+I) - LMIDSZ - 1
      END DO

      CALL MOVEI ( C2PAGE(TRDATC+N+1), RSHIFT-1, C3PAGE(TRDATC+1) )
      CALL MOVEI ( C2PAGE(TRKIDC+N+1), RSHIFT,   C3PAGE(TRKIDC+1) )

C
C     Update the key count in the right child.
C
      C3PAGE ( TRNKC )  =  RSIZE

C
C     The right child is complete.  It's time to update the parent.
C
C     The key at location N in the middle child replaces the left parent
C     key.  The key value is actually re-assigned; the data pointer does
C     move.  The left parent key increases by the number of keys moved
C     into the subtree headed by the left child.
C
      PPAGE ( KEYBAS + LPKIDX )  =  LPKEY + LMIDSZ + 1
      PPAGE ( DATBAS + LPKIDX )  =  C2PAGE( TRDATC + N )

C
C     The parent keys, data pointers, and child pointers at locations
C     LPKIDX+2 onward get shifted left by one position.  The keys lose
C     no  predecessors as the result of these re-arrangements.
C
      DO I = LPKIDX+1, NPKEYS-1
         PPAGE(KEYBAS+I)  =  PPAGE(KEYBAS+I+1)
      END DO

      DO I = LPKIDX+1, NPKEYS-1
         PPAGE(DATBAS+I)  =  PPAGE(DATBAS+I+1)
      END DO

      DO I = LPKIDX+1, NPKEYS
         PPAGE(KIDBAS+I)  =  PPAGE(KIDBAS+I+1)
      END DO

C
C     Zero out the freed locations.
C
      PPAGE ( KEYBAS+NPKEYS   ) = 0
      PPAGE ( DATBAS+NPKEYS   ) = 0
      PPAGE ( KIDBAS+NPKEYS+1 ) = 0

C
C     The only required change to the parent's metadata is
C     updating the key count.  At this point, we can set the
C     underflow flag, depending on the status of the parent.
C
      IF ( PARENT .EQ. ROOT ) THEN

         PPAGE ( TRNKR )  =  PPAGE ( TRNKR ) - 1
         UNDRFL           =  PPAGE ( TRNKR )   .EQ.   0

      ELSE

         PPAGE ( TRNKC )  =  PPAGE ( TRNKC ) - 1
         UNDRFL           =  PPAGE ( TRNKC )   .EQ.   MNKEYC - 1

      END IF

C
C     The last change we must make is to update the node count in
C     the root.
C
      IF ( PARENT .EQ. ROOT ) THEN

         PPAGE ( TRNNOD ) =  PPAGE ( TRNNOD ) - 1

      ELSE
C
C        We won't read in the whole root page; we'll just get the
C        base address of the root and update the affected location.
C
         RBASE  =  ZZEKTRBS ( ROOT )

         CALL DASRDI ( HANDLE, RBASE+TRNNOD, RBASE+TRNNOD, NNODE   )
         CALL DASUDI ( HANDLE, RBASE+TRNNOD, RBASE+TRNNOD, NNODE-1 )

      END IF

C
C     Write out our updates.
C
      CALL ZZEKPGWI ( HANDLE,  PARENT,  PPAGE  )
      CALL ZZEKPGWI ( HANDLE,  LEFT,    C1PAGE )
      CALL ZZEKPGWI ( HANDLE,  RIGHT,   C3PAGE )

C
C     Free the page used by the middle child.
C
      CALL ZZEKPGFR ( HANDLE, INT, MIDDLE )

      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


!     ..........................................................   :
!                                                                  :
!          Block                                                   :
!                       R E C O U P L S                            :
!                                                                  :
!     For Calculation Recoupling Coeficients                       :
!                                                                  :
!     Written by G. Gaigalas,                                      :
!                  Department  of  Computer Science,               :
!                  Vanderbilt University,  Nashville               :
!                                                  February 1994   :
!                                                                  :
!     ..........................................................   :
!
!     --------------------------------------------------------------
!     D L S A 1
!     --------------------------------------------------------------
!                                                                  *
!     Written by G. Gaigalas,                                      *
!     Vanderbilt University,  Nashville            February 1994   *
!
      SUBROUTINE DLSA1(K, JA1, KA, IRE, IAT, REC)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  DOUBLE
      USE CONSTS_C
      USE MEDEFN_C
!...Translated by Pacific-Sierra Research 77to90  4.3E  11:15:35  11/16/01
!...Switches:
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      USE ixjtik_I
      USE sixj_I
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER , INTENT(IN) :: K
      INTEGER , INTENT(IN) :: JA1
      INTEGER  :: KA
      INTEGER , INTENT(IN) :: IRE
      INTEGER , INTENT(OUT) :: IAT
      REAL(DOUBLE) , INTENT(OUT) :: REC
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: IA1, IB1, K1, J1, IT1, IT1S, JJ, IFAZ
      REAL(DOUBLE) :: A1
!-----------------------------------------------
      REC = ZERO
      IA1 = J1QN1(JA1,K) - 1
      IB1 = J1QN2(JA1,K) - 1
      IF (JA1 > 2) THEN
         K1 = IHSH + JA1 - 2
         J1 = J1QN1(K1,K) - 1
         K1 = K1 + 1
         IT1 = J1QN1(K1,K) - 1
         IT1S = J1QN2(K1,K) - 1
      ELSE
         JJ = 1
         IF (JA1 == 1) JJ = 2
         J1 = J1QN1(JJ,K) - 1
         K1 = IHSH + 1
         IT1 = J1QN1(K1,K) - 1
         IT1S = J1QN2(K1,K) - 1
      ENDIF
      IF (IRE == 0) THEN
         IF (IXJTIK(KA,IB1,IA1,J1,IT1,IT1S) == 0) RETURN
         IAT = 1
      ELSE
         CALL SIXJ (KA, IB1, IA1, J1, IT1, IT1S, 0, A1)
         A1 = A1*DSQRT(DBLE((IA1 + 1)*(IT1S + 1)))
         IFAZ = J1 + IT1 + IB1 + KA
         IF ((IFAZ/4)*4 /= IFAZ) A1 = -A1
         REC = A1
         IAT = 1
         IF (JA1 == 1) THEN
            IFAZ = IA1 + IB1 + 2*J1 - IT1 - IT1S
            IF ((IFAZ/4)*4 /= IFAZ) REC = -REC
         ENDIF
      ENDIF
      RETURN
      END SUBROUTINE DLSA1