Program test_euler_sky
INTEGER, PARAMETER :: N_D_FRE=5, N_COEFF=9, DP = KIND (1.d0)
INTEGER :: I_DIAG (N_D_FRE), I_BC (N_D_FRE)
REAL(DP), PARAMETER :: CONST = 96.d0
REAL(DP) :: M_SKY(N_COEFF), S_SKY(N_COEFF), CC_IN (N_D_FRE)
REAL(DP) :: DD (N_D_FRE)
  
  I_DIAG = (/ 1,3,5,7,9 /)
  I_BC   = (/ 1,0,0,0,0 /)
  M_SKY = (/ 2.d0, 1.d0, 4.d0, 1.d0, 4.d0, &
             1.d0, 4.d0, 1.d0, 2.d0 /)
  S_SKY = (/ 1.d0, -1.d0, 2.d0, -1.d0, 2.d0, &
             -1.d0, 2.d0, -1.d0, 1.d0 /)
  S_SKY = CONST * S_SKY
  DD = 0.d0 ; CC_IN = 0.d0

! debug check input
 CALL SKY_FUL_SYM (M_SQ, I_DIAG, CC_IN)
 CALL SKY_FUL_SYM (S_SQ, I_DIAG, DD   )
 
  CALL EULER_SKY_SYM  (N_D_FRE, N_COEFF, I_BC, &
                       I_DIAG, M_SQ, S_SQ, CC_IN, DD)
  print *, 'End test_euler_sky'
contains 

SUBROUTINE EULER_SKY_SYM  (N_D_FRE, N_COEFF, I_BC, &
                           I_DIAG, M_SQ, S_SQ, CC_IN, DD)
 IMPLICIT NONE
 LOGICAL :: FACT, BACK

 INTEGER,  INTENT(IN) :: N_D_FRE, N_COEFF
!            Given Arrays
 INTEGER,  INTENT(IN) :: I_BC   (N_D_FRE)
 INTEGER,  INTENT(IN) :: I_DIAG (N_D_FRE)

 REAL(DP), INTENT(IN) :: CC_IN (N_D_FRE) !b 
 REAL(DP), INTENT(INOUT) :: DD    (N_D_FRE) !b 
 REAL(DP), INTENT(INOUT) :: M_SKY(N_COEFF)             !b
 REAL(DP), INTENT(INOUT) :: S_SKY(N_COEFF)             !b


 INTEGER  :: J, K, T_G, T_S  ! TIME GROUP AND TIME STEP NUMBERS
 INTEGER  :: T_SETS = 1
 REAL(DP) :: START_TIME = 0.d0, TIME_STEP = 0.01d0

! modify for time steping
 S_SQ = M_SQ - TIME_STEP * S_SQ 

 ! apply zero bc
 do 1 k =1, N_D_FRE
   if ( I_BC (k) == 0 ) cycle
     CALL SKY_TYPE_1_SYM (k, 0.d0, S_SQ, CC_IN, I_DIAG)
     CALL SKY_TYPE_1_SYM (k, 0.d0, M_SQ, CC_IN, I_DIAG)
 end do

! factor M_SQ
  FACT = .true. ; BACK = .false.
  call SKY_SOLVE_SYM (M_SQ, CC_IN, X, I_DIAG, FACT, BACK)
! end initial calculations

!     ***  BEGIN TIME INTERGATION GROUPS ***
  DD     = 0.d0 !b
  START_TIME = START_TIME !b
  GROUPS: DO T_G = 1, T_SETS
!    RECOVER TIME INTEGRATION DATA FOR EACH GROUP
     T_STEPS   =  T_STEPS
     TIME_STEP = TIME_STEP
     TS_B      = TS_B

print *, 'beginning euler consistent matrix time integration'
      T_S = 0 ; TIME = START_TIME
      WRITE (N_PRT, 7) T_S, TIME
      7 FORMAT ('STEP = ', I5, 5X, ' TIME = ', 1PE12.5, /, &
                '   DOF           VALUE')
      WRITE (N_PRT, 8) (K, DD (K), K = 1, N_D_FRE)
      8 FORMAT (I6, 2X, 1PE12.5)


!         *** TIME STEPPING WITHIN EACH TIME GROUP ***
     STEPS: DO T_S = 1, T_STEPS
       TIME = START_TIME + T_S * TIME_STEP
       WRITE (N_PRT, 7) T_S, TIME

!b call forcer
      CALL VECTOR_MULT_SKY_SYM (S_SQ, DD, I_DIAG, RESULT, N_D_FRE, N_COEFF)
      CC = CC_IN *  TIME_STEP + RESULT
 ! apply zero bc
      do 1 k =1, N_D_FRE
        if ( I_BC (k) == 0 ) cycle
          CC (k) = 0.d0
      end do

! solve for solution at time t
  FACT = .false. ; BACK = .true.
  call SKY_SOLVE_SYM (M_SQ, CC, DD, I_DIAG, FACT, BACK)

!-->      *** SOLVE FOR UNKNOWN NODAL PARAMETERS ***
      FACT = .FALSE. ; BACK = .TRUE.
      CALL  SKY_SOLVE_SYM (M_SKY, CC, DD, I_DIAG, FACT, BACK)
      WRITE (N_PRT, 8) (K, DD (K), K = 1, N_D_FRE)

!          *** PRINT RESULTS ***
     END DO STEPS
   END DO GROUPS
END SUBROUTINE EULER_SKY_SYM

SUBROUTINE SKY_TYPE_1_SYM (LOC, VALUE, S, C, I_DIAG, &  !b
                           N_D_FRE, N_COEFF)  !b
!b SUBROUTINE SKY_TYPE_1_SYM (LOC, VALUE, S, C, I_DIAG)
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
! APPLY TYPE 1 MODIFICATION TO SYMMERTIC SKYLINE EQS
!      S*D = C,  D(LOC) = VALUE
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
!b Use System_Constants
 IMPLICIT NONE
 INTEGER,  INTENT(IN)    :: N_D_FRE, N_COEFF  !b
 INTEGER,  INTENT(IN)    :: I_DIAG (N_D_FRE), LOC
 REAL(DP), INTENT(IN)    :: VALUE
 REAL(DP), INTENT(INOUT) :: S (N_COEFF), C (N_D_FRE)
 INTEGER   :: I, ID, INV, ITOP

! N_D_FRE = NUMBER OF EQUATIONS
! N_COEFF = NUMBER OF NON-ZERO COEFFICIENTS IN S
! LOC     = DOF NUMBER OF CONSTRAINED PARAMETER
! VALUE   = GIVEN VALUE OF DOF NUMBER LOC
! S       = SYSTEM SQUARE MATRIX IN SKYLINE STORAGE MODE
! C       = FULL COLUMN MATRIX
! I_DIAG  = POINTER TO DIAGONAL COEFFICIENT IN S
!
  IF ( LOC < 1 .OR. LOC > N_D_FRE ) STOP 'BAD LOC, SKY_TYPE_1_SYM'
!  SUBTRACT COLUMN*VALUE FROM RHS
  DO I = 1, N_D_FRE
!    FIND S(I, LOC) IN S VECTOR
    ID = MAX0 (I, LOC)
    INV = I_DIAG (ID) - IABS (I - LOC)
    ITOP = 1
    IF (ID > 1) ITOP = I_DIAG (ID-1) + 1
!    IS IT OUTSIDE SKYLINE AND THUS ZERO?
    IF (INV >= ITOP) THEN
      C (I) = C (I) - VALUE * S (INV)
      S (INV) = 0.d0
    ENDIF
  END DO
!  RESET THE EQUATION ROW
  S (I_DIAG (LOC)) = 1.d0
  C (LOC) = VALUE
END SUBROUTINE SKY_TYPE_1_SYM

SUBROUTINE VECTOR_MULT_SKY_SYM (S_SKY, VECTOR, I_DIAG, RESULT, &
                                N_D_FRE, N_COEFF)
!b SUBROUTINE VECTOR_MULT_SKY_SYM (S_SKY, VECTOR, I_DIAG, RESULT)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! MULTIPLY A SYMMETRIC SKYLINE ARRAY BY A VECTOR, SKY*VECTOR = RESULT
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!b Use System_Constants ! N_D_FRE, N_COEFF, NOT_SYM, MAX_NP, SYMMETRIC
!b Use Interface_Header
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: N_D_FRE, N_COEFF !b
 REAL(DP), INTENT(IN)  :: S_SKY       (N_COEFF)     ! always used
 REAL(DP), INTENT(IN)  :: VECTOR      (N_D_FRE)     ! always used
 INTEGER,  INTENT(IN)  :: I_DIAG      (N_D_FRE)     ! always used
 REAL(DP), INTENT(OUT) :: RESULT      (N_D_FRE)     ! always used
!          Automatic Arrays
 INTEGER   :: ROW              ! equation number of dof at node
 INTEGER   :: IN_SKY           ! location in skyline, 0 if absent
 INTEGER   :: COL, COL_START, COL_STOP   ! column search at row
 REAL (DP) :: DOT                        ! row dot product

! I_DIAG  = DIAGONAL LOCATION IN SKYLINE VECTOR
! N_COEFF = NUMBER OF TERMS IN UPPER STIFFNESS
! N_G_DOF = NUMBER OF DOF PER NODE
! N_D_FRE = TOTAL NUMBER OF EQUATIONS
! NOT_SYM = NUMBER OF TERMS IN LOWER STIFFNESS, 0 OR N_COEFF
! S_SKY   = UPPER HALF OF SS MATRIX IN SKYLINE MODE

! NOTE: Alternate way is to loop over non_zeros in skyline and
! multiply by the corresponding vector

  DO ROW = 1, N_D_FRE                 !  LOOP OVER EQUATIONS
!       FIND COLUMN NUMBER BOUNDS FOR NON-ZERO COEFFICIENTS IN ROW
    COL_START = FIRST_COL_AT_SKY_ROW (ROW, I_DIAG, N_D_FRE) !b
    COL_STOP  =  LAST_COL_AT_SKY_ROW (ROW, I_DIAG, N_D_FRE) !b
!     CORRECT FOR ANY INITIALLY NULL ROWS (GRAB DIAGONAL)
    IF ( COL_START > ROW ) COL_START = ROW
    IF ( COL_STOP  < ROW ) COL_STOP  = ROW
!     MULTIPLY ROW OF NON_ZEROS IN SYSTEM EQUATIONS
    DOT = 0.d0
    DO COL = COL_START, COL_STOP                   ! column loop
      IN_SKY = LOCATION_IN_SKY (ROW, COL, I_DIAG, N_D_FRE)  ! locate
      IF ( IN_SKY > 0 ) DOT = DOT + S_SKY (IN_SKY) * VECTOR (COL)
    END DO ! OVER COLUMNS
    RESULT (ROW) = DOT
  END DO ! OVER ROWS
END SUBROUTINE VECTOR_MULT_SKY_SYM

FUNCTION FIRST_COL_AT_SKY_ROW (ROW, I_DIAG, N_D_FRE) RESULT (COLUMN)
!b FUNCTION FIRST_COL_AT_SKY_ROW (ROW, I_DIAG) RESULT (COLUMN)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!     FIND FIRST COLUMN IN A ROW IN A SKYLINE STORAGE
!           WARNING: CAN OCCUR AFTER A DIAGONAL
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!   (ALSO FIRST ROW IN A COLUMN IN A SKYLINE STORAGE)
!b Use system_constants ! N_D_FRE
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: N_D_FRE
 INTEGER, INTENT(IN) :: ROW               ! row (or column) number
 INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE)  ! skyline diagonals
 INTEGER :: COLUMN          ! first column ( or row)
 INTEGER :: COL             ! loop
 INTEGER :: LOCATION_IN_SKY

! N_D_FRE = TOTAL NUMBER OF EQUATIONS
    COLUMN = 0                           ! may be null
    DO COL = 1, N_D_FRE                  ! forward column loop
      IF ( LOCATION_IN_SKY (ROW, COL, I_DIAG) > 0 ) THEN ! non-zero
        COLUMN = COL ; EXIT              ! search over
      END IF
    END DO
    IF ( COLUMN == 0 ) THEN
      PRINT *,'WARNING, NULL FIRST_COL_AT_SKY_ROW'
      N_WARN = N_WARN + 1 ! INCREMENT WARNING
    END IF
END FUNCTION FIRST_COL_AT_SKY_ROW

FUNCTION LAST_COL_AT_SKY_ROW (ROW, I_DIAG, N_D_FRE) RESULT (COLUMN)
!b FUNCTION LAST_COL_AT_SKY_ROW (ROW, I_DIAG) RESULT (COLUMN)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!    FIND LAST COLUMN IN A ROW IN A SKYLINE STORAGE
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!   (ALSO LAST ROW IN A COLUMN IN A SKYLINE STORAGE)
!b Use system_constants ! N_D_FRE
! Use Interface_Header
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: N_D_FRE
 INTEGER, INTENT(IN) :: ROW               ! row (or column) number
 INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE)  ! skyline diagonals
 INTEGER :: COLUMN          ! first column ( or row)
 INTEGER :: COL             ! loop
 INTEGER :: LOCATION_IN_SKY

! N_D_FRE = TOTAL NUMBER OF EQUATIONS
    COLUMN = 0                        ! maybe a null row
    DO COL = N_D_FRE, 1, -1           ! backward column loop
      IF ( LOCATION_IN_SKY (ROW, COL, I_DIAG) > 0 ) THEN ! non-zero
        COLUMN = COL ; EXIT           ! search over
      END IF
    END DO
    IF ( COLUMN == 0 ) THEN
      PRINT *,'WARNING, NULL LAST_COL_AT_SKY_ROW'
      N_WARN = N_WARN + 1 ! INCREMENT WARNING
    END IF
END FUNCTION LAST_COL_AT_SKY_ROW

FUNCTION LOCATION_IN_SKY (I, J, I_DIAG, N_D_FRE) RESULT (IN_SKY)
!b FUNCTION LOCATION_IN_SKY (I, J, I_DIAG) RESULT (IN_SKY)
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
! CONVERT (I,J) FULL SYMMETRIC MATRIX SUBSCRIPTS TO
! (IN_SKY) SUBSCRIPT OF VECTOR SKYLINE STORAGE MODE.
! IN_SKY = 0 IF OUTSIDE SKYLINE AND THUS HAS 0 VALUE.
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
! ASSUMING COLUMNS STORED FROM TOP DOWN
!b Use system_constants ! N_D_FRE
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: N_D_FRE
 INTEGER, INTENT(IN)  :: I, J            ! full subscripts
 INTEGER, INTENT(IN)  :: I_DIAG(N_D_FRE)  ! diagonal locations
 INTEGER              :: IN_SKY          ! location in skyline
 INTEGER              :: ID              ! max equation

! N_D_FRE   = TOTAL NO OF SYSTEM EQUATIONS
    ID  = MAX0 (I,J)                      ! max equation
    IF ( ID < 1 .OR. ID > N_D_FRE) STOP 'BAD ID, LOCATION_IN_SKY'
    IN_SKY = I_DIAG(ID) - IABS(I-J)        ! possible non-zero
!     TEST FOR ZERO ENTRY OUTSIDE SKYLINE
    IF ( ID > 1 ) THEN
      IF ( IN_SKY <= I_DIAG(ID - 1) ) IN_SKY = 0
    END IF
END FUNCTION LOCATION_IN_SKY

Program test_euler_sky
