! copyright 2005,  J. E. Akin, all rights reserved.
SUBROUTINE ASSEMBLE_SKYLINE_EQS (I_DIAG, NODES, S_SKY, CC, X, DD_OLD, &
                                 S_SKY_LOWER, ITER, M_SKY, M_SKY_LOWER)
                              !b S_SKY_LOWER, ITER, M_SKY)
                          !b X, DD_OLD, S_SKY_LOWER, ITER, I_BC) !b
                          !b X, DD_OLD, S_SKY_LOWER, ITER)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!           ASSEMBLE SYSTEM EQUATIONS AND STORE POST
!                     SOLUTION ELEMENT DATA
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data       ! for PT (LT_PARM, LT_QP), WT (LT_QP),
                         ! G (LT_GEOM, LT_QP), H (LT_FREE),
                         ! DLG (LT_PARM, LT_GEOM, LT_QP),
                         ! DLH (LT_PARM, LT_FREE, LT_QP),
                         ! S (LT_FREE, LT_FREE), C (LT_FREE),
                         ! EL_M (LT_FREE, LT_FREE), DIAG_M (LT_FREE),
                         ! D (LT_FREE), ELEM_NODES (LT_N)
Use Sys_Properties_Data
Use Interface_Header
 IMPLICIT NONE
!                   ALWAYS USED
 INTEGER,  INTENT(IN)    :: ITER  ! Current iteration number
 INTEGER,  INTENT(IN)    :: NODES (L_S_TOT, NOD_PER_EL)
 INTEGER,  INTENT(IN)    :: I_DIAG (N_D_FRE)
 !b INTEGER,  INTENT(IN)    :: I_BC   (MAX_NP)     !b
 REAL(DP), INTENT(IN)    :: X (MAX_NP, N_SPACE) !b
 REAL(DP), INTENT(IN)    :: DD_OLD (N_D_FRE)    !b
 REAL(DP), INTENT(INOUT) :: CC (N_D_FRE)
 REAL(DP), INTENT(INOUT) :: S_SKY (N_COEFF)
 REAL(DP), INTENT(INOUT) :: S_SKY_LOWER (NOT_SYM)
 REAL(DP), INTENT(INOUT) :: M_SKY (M_COEFF)
 REAL(DP), INTENT(INOUT) :: M_SKY_LOWER (NOT_SYM) ! effective unsym mass
!                 Locals
 INTEGER  :: IE, J, LT
 REAL(DP) :: COL_SUMS   ! Error checking

!   Work items for time integration
 REAL(DP)  :: K (LT_FREE, LT_FREE), M (LT_FREE, LT_FREE)
 REAL(DP)  :: F (LT_FREE)

!                   VARIABLES:
! CC         = SYSTEM EQUATIONS COLUMN MATRIX
! COL_SUMS   + SUM OF ALL COLUMN MATRICES
! D          = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT
! DD_OLD     = SYSTEM NODAL PARAMETERS FROM LAST ITERATION
! DELETED    = TRUE IF THE ELEMENT HAS BEEN DELETED
! DGH        = GLOBAL DERIVATIVES INTERPOLATION FUNCTIONS
! DIAG_M     = DIAGONALIZED ELEMENT MASS MATRIX (LT_FREE)
! DLG        = LOCAL DERIVATIVES GEOMETRIC INTERPOLATION
! DLH        = LOCAL DERIVATIVES INTERPOLATION FUNCTIONS
! E          = CONSTITUTIVE MATRIX
! EB         = PRODUCT OF E*B
! ELEM_NODES = THE NOD_PER_EL INCIDENCES OF THE ELEMENT
! EL_M       = ELEMENT MASS MATRIX
! FLT_EL     = REAL PROPERTIES OF ELEMENTS
! FLT_MISC   = MISCELLANEOUS REAL PROPERTIES OF SYSTEM
! FLT_NP     = REAL PROPERTIES OF SYSTEM NODES
! G          = GEOMETRIC INTERPOLATION FUNCTIONS
! H          = SOLUTION INTERPOLATION FUNCTIONS
! H_INTG     = INTEGRAL OF INTERPOLATION FUNCTIONS
! ITER       = CURRENT ITERATION NUMBER (not used)
!b ! I_BC       = NODAL BOUNDARY RESTRAINT INDICATOR ARRAY  !b
! I_DIAG     = DIAGONAL LOCATION IN SKYLINE VECTOR
! INDEX      = SYSTEM DOF NUMBERS ASSOCIATED WITH ELEMENT
! L_HOMO     = 1, IF ELEMENT PROPERTIES ARE HOMOGENEOUS
! L_S_TOT    = TOTAL NUMBER OF ELEMENTS & THEIR SEGMENTS
! LP_FIX     = SYSTEM ARRAY OF INTEGER ELEM PROPERTIES
! L_PROP     = ARRAY INTEGER ELEMENT PROPERTIES
! LP_TEST    > 0, IF ELEMENT PROPERTIES HAVE BEEN DEFINED
! L_SHAPE    = SHAPE FLAG FOR QUADRATURE RULE SELECTION
! LT_DATA    = ARRAY OF CONSTANTS FOR EACH ELEMENT TYPE
! LT         = ELEMENT TYPE NUMBER (IF USED)
! LT_QP      = NUMBER OF QUADRATURE PTS FOR ELEMENT TYPE
! LT_USER    = OPTIONAL USER CODE ASSOCIATED WITH TYPE
! L_TYPE     = ELEMENT TYPE NUMBER
! MAX_NP     = NUMBER OF SYSTEM NODES
! MODE       = STORAGE FLAG, 0=SYMMETRIC, OR 1=UNSYMMETRIC SKYLINE
! MISC_FIX   = MISCELLANEOUS INTEGER SYSTEM PROPERTIES
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! N_COEFF    = TOTAL NUMBER OF TERMS IN SS
! N_D_FRE    = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM
! N_ELEMS    = NUMBER OF ELEMENTS
! N_EL_FRE   = NUMBER OF DEGREES OF FREEDOM PER ELEMENT
! N_G_DOF    = NUMBER OF PARAMETERS PER NODE
! N_GEOM     = NUMBER OF GEOMETRY NODES
! N_HOMO     = 1, IF NODAL PROPERTIES ARE HOMOGENEOUS
! N_ITER     = NO. OF ITERATIONS TO BE RUN (USUALLY 1)
! N_MAT      = NUMBER OF MATERIAL TYPES
! NODES      = NODAL INCIDENCES OF ALL ELEMENTS
! N_PARM     = DIMENSION OF PARAMETRIC SPACE
! NP_FIX     = INTEGER PROPERTIES AT ALL NODES
! N_QP       = NUMBER OF QUADRATURE POINTS, >= LT_QP
! N_R_B      = NUMBER OF ROWS IN B AND E MATRICES
! N_SPACE    = DIMENSION OF SPACE
! N_FILE1    = UNIT FOR POST SOLUTION MATRICES STORAGE
! N_FILE2    = OPTIONAL UNITS FOR USER (USED WHEN > 0)
! N_FILE3,4  = OPTIONAL UNITS FOR USER (USED WHEN > 0)
! NUL_COL    = 1, IF ELEMENT COLUMN MATRIX IS ALWAYS ZERO
! PT         = QUADRATURE COORDINATES
! S          = ELEMENT SQUARE MATRIX
! S_SKY      = SYSTEM EQUATIONS SQUARE MATRIX
! STRAIN     = STRAIN OR GRADIENT VECTOR
! STRAIN_0   = INITIAL STRAIN OR GRADIENT VECTOR
! STRESS     = STRESS VECTOR
! WT         = QUADRATURE WEIGHTS
! X          = COORDINATES OF SYSTEM NODES
! XYZ        = SPACE COORDINATES AT A POINT
 
! GENERATE ELEMENT EQUATIONS & POST SOLUTION MATRICES
  WRITE (N_PRT, "(/, 'BEGINNING ASSEMBLY')")
  IF ( I_BUG > 0 .AND. N_ITER > 1 ) PRINT *, 'ITERATION ', ITER
  IF ( U_FLUX  > 0 ) REWIND U_FLUX  
  IF ( N_FILE1 > 0 ) REWIND N_FILE1 
  IF ( N_FILE2 > 0 ) REWIND N_FILE2 
  COL_SUMS = 0.d0

  LT = 1 ; LAST_LT = 0 ! INITIALIZE ELEMENT TYPE

!     LOOP OVER ALL STANDARD ELEMENTS FIRST, 
!     THEN ANY OF THEIR BOUNDARY SEGMENT ELEMENTS

  WRITE (N_PRT, "('BEGINNING STANDARD ELEMENT ASSEMBLY')")
  DO IE = 1, L_S_TOT ! for elements and any boundary segments
!b call at(125)
    IF ( DELETED (IE) ) CYCLE ! to skip deleted element
    CALL SET_THIS_ELEMENT_NUMBER (IE)
    IF ( N_MIXED > 0 .AND. IE == (N_ELEMS + 1) ) WRITE (N_PRT, &
      "('BEGINNING MIXED_BC SEGMENTS ASSEMBLY')")
    IF ( N_F_SEG > 0 .AND. IE == (N_ELEMS + N_MIXED + 1) ) &
      WRITE (N_PRT, "('BEGINNING FLUX BOUNDARY SEGMENTS ASSEMBLY')")

!-->     GET ELEMENT TYPE NUMBER
    IF ( N_L_TYPE > 1 ) LT = L_TYPE (IE)
!    SAME AS LAST TYPE ?
    IF ( LT /= LAST_LT ) THEN ! this is a new type, or first element
!b call at(137)
!      GET CONTROLS FOR THIS TYPE, AND ALLOCATE ITS ARRAYS
      CALL SET_ELEM_TYPE_INFO (LT) ! including LAST_LT
    END IF ! a new element type

!        INITIALIZE AND CHECK IF ACTIVE
    !b    call list_type_alloc_status !b
    S = 0.d0 ; EL_M = 0.d0 ; C = 0.d0

!-->     EXTRACT ELEMENT NODE NUMBERS
      ELEM_NODES = GET_ELEM_NODES (IE, LT_N, NODES)    

!-->     CALCULATE DEGREE OF FREEDOM NUMBERS
      INDEX = GET_ELEM_INDEX (LT_N, ELEM_NODES)   

if (debug_el_sq .AND. (N_ITER > 0 .OR. DYNAMIC .OR. TRANSIENT)) then
   !b call at(155) ; 
   print *,'DD_OLD ' ; call rprint(DD_OLD, 1, N_D_FRE, 1)
end if
!-->     GENERATE ELEMENT PROBLEM DEPENDENT MATRICES
      IF ( EXAMPLE > 0 ) THEN ! Select an old example
        CALL GENERATE_EL_MATRICES_NEW (IE, X, DD_OLD)
        IF ( DEBUG_EL_COL .OR. DEBUG_EL_SQ ) THEN
          PRINT *, 'Exiting GENERATE_EL_MATRICES_NEW eith C, S:'
          IF ( IE <= 3 )  CALL RPRINT (C, 1, LT_FREE, 1) !b
          IF ( IE <= 3 )  CALL RPRINT (S, LT_FREE, LT_FREE, 1) !b
        END IF
      ELSE
        CALL GENERATE_EL_MATRICES     (IE, X, DD_OLD)
      END IF ! Select an old example
!b call at(153)

!-->     SAVE ELEMENT MATRICES FOR ELEMENT REACTIONS
      IF ( L_REACT > 0 ) WRITE (L_REACT) S, C

!-->     STORE THE MATRICES IN SYSTEM EQUATIONS
      IF ( STATIC ) THEN
        IF ( I_BUG > 0 ) PRINT *,'BEGIN STATIC ANALYSIS'
        CALL  SKY_STORE_SQ (LT_FREE, INDEX, I_DIAG, S, S_SKY, &
                            S_SKY_LOWER)
        IF ( NUL_COL == 0) THEN ! SOME NON_ZERO ELEMENT SOURCES
          CALL STORE_COLUMN (LT_FREE, INDEX, C, CC)
          COL_SUMS = COL_SUMS + SUM ( C )  ! Error checking
          IF ( DEBUG_EL_COL .AND. IE <= 3 )  CALL RPRINT (C, 1, LT_FREE, 1) !b
        ELSE ! validate key words
          IF ( IE == 1 .AND. ANY (C /= 0.d0 ) ) THEN
            N_WARN = N_WARN + 1
            PRINT *, 'WARNING: KEYWORD el_no_col IS WRONG'
          END IF
        END IF ! ELEMENT SOURCES

      ELSE IF ( TRANSIENT ) THEN ! BRANCH ON TIME INTEGRATION METHOD 

!       So element matrix     S * D (t) + EL_M * D_dot (t) = C (t)
!       Converts to algebraic K * D (t) = F (t) with
!       K a function of S, EL_M, STEP_SIZE
!       M a function of S, EL_M, STEP_SIZE, so that
!       F is a function of M * D and C

        IF ( I_BUG > 0 ) PRINT *,'BEGIN TRANSIENT ANALYSIS'

!       Convert element level arrays
!               S (LT_FREE, LT_FREE), EL_M (LT_FREE, LT_FREE),
!               C (LT_FREE), or DIAG_M  (LT_FREE) with
!               D (LT_FREE) to form
!       One step transient forms
!               K (LT_FREE, LT_FREE), M (LT_FREE, LT_FREE),
!               F (LT_FREE) for assembly into 
!               K(S, EL_M, Dt) * D(t) = F(C, S, EL_M, Dt, D)
        CALL ONE_STEP_TRANSIENT_AT_EL (K, M, F)

!           ASSEMBLE, ASSUMING CONSTANT C W.R.T. TIME
        CALL  SKY_STORE_SQ (LT_FREE, INDEX, I_DIAG, K, S_SKY, & !b
                            S_SKY_LOWER)
        CALL  SKY_STORE_SQ (LT_FREE, INDEX, I_DIAG, M, M_SKY, & !b
                            M_SKY_LOWER)
        CALL STORE_COLUMN (LT_FREE, INDEX, F, CC) !b

!         ASSEMBLE REACTION ROWS FILE IF TRANSIENT REACTIONS DESIRED
!b         like SAVE_SKY_ROWS_GEN for use in GET_REACTIONS

        IF ( DEBUG_ASSEMBLY ) THEN
          print *,'K_SKY at 214 for element ', THIS_EL
          CALL SKY_FUL_SYM (S_SKY, I_DIAG, CC)
          print *,'M_SKY at 216'
          CALL SKY_FUL_SYM (M_SKY, I_DIAG, CC)
        END IF
        IF ( DEBUG_EL_COL ) THEN
          PRINT *,'EL, F, CC ', THIS_EL
          CALL RPRINT (F, 1, LT_FREE, 1) !b
          CALL RPRINT (CC, 1, N_D_FRE, 1) !b
        END IF

      ELSE IF ( DYNAMIC ) THEN
        IF ( I_BUG > 0 ) PRINT *,'BEGIN DYNAMIC ANALYSIS'

!           ASSEMBLE, ASSUMING CONSTANT C W.R.T. TIME
        CALL  SKY_STORE_SQ (LT_FREE, INDEX, I_DIAG, S, S_SKY, & !b
                            S_SKY_LOWER)
!          Symmetric mass matrix only
        CALL  SKY_STORE_SQ_SYM (LT_FREE, INDEX, I_DIAG, EL_M, M_SKY)
!          Symmetric element damping matrix only
        IF ( EL_DAMPING ) CALL  SKY_STORE_SQ_SYM (LT_FREE, INDEX, &
                                  I_DIAG, DAMP, M_SKY_LOWER)
        CALL STORE_COLUMN (LT_FREE, INDEX, C, CC) !b

      ELSE IF ( EIGEN ) THEN
          !b PRINT *, 'BEGINNING EIGEN ANALYSIS'
!           ASSEMBLE STIFFNESS AND MASS MATRICES FROM ELEMENT
        !b CALL  SKY_STORE_SQ (LT_FREE, INDEX, I_DIAG, S, S_SKY, &
                            !b S_SKY_LOWER)
        !b CALL  SKY_STORE_SQ (LT_FREE, INDEX, I_DIAG, EL_M, M_SKY, &
                            !b M_SKY_LOWER)
          STOP 'ASSEMBLE_SKYLINE_EQS: EIGEN ASSEMBLY UNDEFINED'

      ELSE ! UNKNOWN DEFAULT
          STOP 'ASSEMBLE_SKYLINE_EQS: ELEMENT ASSEMBLY CLASS UNDEFINED'
      END IF ! ANALYSIS CLASS

      !b IF ( .NOT. STATIC ) THEN ! NEED MASS ASSEMBLY
        !b IF ( DIAGONAL_MASS ) THEN
          !b CALL  DIAGONALIZE_SQ_MATRIX (LT_FREE, EL_M, DIAG_M)
          !b CALL  STORE_COLUMN (LT_FREE, INDEX, DIAG_M, M_SKY)
        !b ELSE
          !b CALL  SKY_STORE_SQ_SYM (LT_FREE, INDEX, I_DIAG, EL_M, M_SKY)
          !b IF ( .NOT. SYMMETRIC ) THEN
            !b WRITE (N_PRT, *) 'WARNING: '
            !b N_WARN = N_WARN + 1
          !b END IF
        !b END IF ! DIAGONAL
      !b END IF ! NOT STATIC

      IF ( IE == N_ELEMS ) WRITE (N_PRT,                 &
                "('COMPLETED STANDARD ELEMENT ASSEMBLY')")
      IF ( N_MIXED > 0 .AND. IE == (N_ELEMS + N_MIXED) )        &
        WRITE (N_PRT, "('COMPLETED MIXED_BC SEGMENTS ASSEMBLY')")
      IF ( N_F_SEG > 0 .AND.                                   &
           IE == (N_ELEMS + N_MIXED + N_F_SEG) ) WRITE (N_PRT, &
           "('COMPLETED FLUX BOUNDARY SEGMENTS ASSEMBLY')")
      IF ( IE == L_S_TOT .AND. (N_MIXED > 0 .OR. N_F_SEG > 0) ) &
        WRITE (N_PRT, "('ASSEMBLY COMPLETE.', /)")
  END DO ! over all elements and boundary segments
  IF ( NUL_COL == 0 .AND. COL_SUMS == 0.d0 ) THEN
    PRINT *, 'WARNING: ALL ELEMENT SOURCE VECTORS ARE NULL,'
    PRINT *, 'PUT el_no_col IN KEYWORD CONTROL'            
    N_WARN = N_WARN + 1 ! INCREMENT WARNING
  END IF

  IF ( DEBUG_EL_COL .AND. IE <= 3 )  CALL RPRINT (CC, 1, N_D_FRE, 1) !b
  WRITE (N_PRT, "(/, 'ASSEMBLY COMPLETED')")
END SUBROUTINE ASSEMBLE_SKYLINE_EQS
                               
SUBROUTINE CEQ_SKYLINE (N_RES_EQ, NDEX_C_EQ, I_DOF_HI, I_DIAG)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! FIND SKYLINE STORAGE REQUIRED BY CONSTRAINT EQUATION
! MODIFICATION PROCEDURES.  IT ALSO UPDATES I_DIAG AND N_COEFF
! AND I_DOF_HI.  SOME OF THE CALCULATIONS ARE CONSERVATIVE
! (MORE STORAGE PROVIDED THAN NEEDED)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Interface_Header
 IMPLICIT NONE
 INTEGER, INTENT(INOUT) :: I_DOF_HI (N_D_FRE), I_DIAG (N_D_FRE)
 INTEGER, INTENT(IN)    :: NDEX_C_EQ (MAX_ACT, N_CEQ)
 INTEGER, INTENT(IN)    :: N_RES_EQ (MAX_ACT)
 INTEGER :: IC, IEQ, II, IIC, IICTOP, IJ, IJTOP, JJ
 INTEGER :: MINDEP, NR,  N_R_BOT, NRTOP, NTEST, ILOOP

! MAX_ACT    = NO ACTIVE CONSTR TYPES
! N_CEQ      = TOTAL NO CONSTR EQS
! NDEX_C_EQ  = (I,J) CONSTR DOF NO I OF EQ J
! N_D_FRE    = TOTAL NO OF SYSTEM DEGREES OF FREEDOM

  IF ( I_BUG == 1 ) THEN
    WRITE (N_BUG, * ) "SKYLINE, I_DOF_HI: ", &
      (I_DOF_HI (II),  II = 1, N_D_FRE)
    WRITE (N_BUG, * ) "I_DIAG: ", &
      (I_DIAG (II),  II = 1, N_D_FRE)
  END IF
  IEQ = N_RES_EQ (1)
!-->  LOOP OVER NON DIAGONAL CONSTRAINTS
  DO IC = 2, MAX_ACT
    NTEST = N_RES_EQ (IC)
    IF ( NTEST > 0 ) THEN
!-->     * LOOP OVER TYPE IC EQUATIONS
      DO ILOOP = 1, NTEST
        IEQ = IEQ + 1
!
!       * THE MODIFICATIONS OF STORAGE ARE DONE FOR THE CURRENT
!       * EQUATION: DJ + SUM ( AK * DK ) = E
!       * TWO STEPS ARE NEEDED FOR THE STORAGE MODIFICATION
!
!       * FIND THE NR DOF AND LOWEST "MINDEP" OF THE DEPENDENT DOFS
        NR = NDEX_C_EQ (1, IEQ)
        MINDEP = N_D_FRE
        DO JJ = 2, IC
          MINDEP = MIN (NDEX_C_EQ (JJ, IEQ), MINDEP)
        END DO
!
!       * FIND THE TOP AND THE BOTTOM OF THE SKYLINE OF COLUMN NR 
        CALL SKY_TOP_BOT (I_DIAG, NR, NRTOP, N_R_BOT)
!
!       * STEP 1: HORIZONTAL FILL-IN (READ THEORY)
!
        DO IJ = NRTOP, N_R_BOT
          IJTOP = IJ - I_DOF_HI (IJ) + 1
!          * SEARCH FOR "PROBABLE NON ZEROS" IN COL NR & MODIFY TOP
          IF ( IJTOP <= NR .OR. IJ < NR ) THEN
            IJTOP = MIN (MINDEP, IJTOP)
            I_DOF_HI (IJ) = IJ - IJTOP + 1
          END IF
        END DO
!
!       * STEP 2: MODIFY THE CONSTRAINT DOFS COLS
!       * THE TOP WILL BE THE MIN OF THREE NUMBERS
!
        DO II = 1, IC
          IIC    = NDEX_C_EQ (II, IEQ)
          IICTOP = IIC - I_DOF_HI (IIC) + 1
          IICTOP = MIN (IICTOP, MINDEP)
          IICTOP = MIN (IICTOP, NRTOP)
          I_DOF_HI (IIC) = IIC - IICTOP + 1
        END DO
!
!       * END OF MODIFICATIONS OF THE SKYLINE STORAGE
!       * UPDATE I_DIAG AND N_COEFF
!
        I_DIAG  = GET_SKY_DIAG (I_DOF_HI) 
        N_COEFF = I_DIAG (N_D_FRE)
!
        IF ( I_BUG == 1 ) THEN
          WRITE (N_BUG, * ) "SKYLINE, I_DOF_HI: ", &
            (I_DOF_HI (II),  II = 1, N_D_FRE)
          WRITE (N_BUG, * ) "I_DIAG: ", &
            (I_DIAG (II),  II = 1, N_D_FRE)
        END IF
      END DO
    END IF
  END DO
END SUBROUTINE CEQ_SKYLINE

SUBROUTINE MESH_SKYLINE (NODES, I_DOF_HI, I_DIAG)
! * * * * * * * * * * * * * * * * * * * * * * * * *
! USE MESH COLUMN HEIGHTS TO FIND DIAGONAL COEFFICIENTS
! FOR SYMMETRIC SKYLINE STORAGE MODE, ASSUMING SYMMETRIC
! COLUMNS STORED FROM TOP DOWN
! * * * * * * * * * * * * * * * * * * * * * * * * *
!  TOTAL NUMBER OF SQ MATRIX TERMS = I_DIAG(N_D_FRE)
Use System_Constants
Use Interface_Header
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: NODES (L_S_TOT, NOD_PER_EL)
 INTEGER, INTENT(OUT) :: I_DIAG (N_D_FRE), I_DOF_HI (N_D_FRE)
 INTEGER :: J, COUNT, N_PT   ! LOOP, INACTIVE COUNTER

! I_DOF_HI(I)  = COLUMN HEIGHT OF EQ I, WITH DIAG
! I_DIAG(I)   = LOCATION OF DIAG OF I-TH EQ NUMBER
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! N_D_FRE    = TOTAL NO OF SYSTEM EQUATIONS
! N_ELEMS    = NUMBER OF ELEMENTS
! NODES      = NODAL INCIDENCES OF ALL ELEMENTS

  CALL SKY_HI (NODES, I_DOF_HI)
  !b IF ( I_BUG > 0 ) PRINT *, 'I_DOF_HI, ', I_DOF_HI 
!  TEST FOR INACTIVE EQUATIONS TO BE CORRECTED
  COUNT = 0
  DO J = 1, N_D_FRE
    IF ( I_DOF_HI (J) <  0 ) STOP 'NEGATIVE HEIGHTS, MESH_SKYLINE'
    IF ( I_DOF_HI (J) == 0 ) THEN ! INACTIVE EQ WITH NO ENTRIES
      IF ( COUNT == 0 ) THEN
        PRINT *,'WARNING: NULL EQUATION AT EQUATION NUMBER = ', J  
        CALL GET_NODE_AT_DOF (J, N_G_DOF, N_PT)
        PRINT *,'NOTE: FIRST NULL EQUATION IS AT NODE ', N_PT
        N_WARN = N_WARN + 1 ! INCREMENT WARNING
      END IF
      COUNT        = COUNT + 1
      I_DOF_HI (J) = 1  
    END IF  
  END DO 
  IF ( COUNT > 0 ) THEN
    PRINT *, 'WARNING: TOTAL NUMBER OF NULL EQS = ', COUNT
    N_WARN = N_WARN + 1 ! INCREMENT WARNING
  END IF
!  CONVERT TO STORAGE LOCATIONS
  I_DIAG  = GET_SKY_DIAG (I_DOF_HI) 
END SUBROUTINE MESH_SKYLINE
                           
SUBROUTINE SKY_BC_GEN (N_RES_EQ, COEF_C_EQ, NDEX_C_EQ, S_SKY, &
                       CC, I_DIAG, S_SKY_LOWER)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!   APPLY BOUNDARY CONSTRAINT EQUATIONS TO GENERAL SKYLINE SYSTEM
!                         OF EQS, S*D = CC
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 REAL(DP), INTENT(IN)    :: COEF_C_EQ (MAX_ACT, N_CEQ)
 INTEGER,  INTENT(IN)    :: N_RES_EQ  (MAX_ACT),             &
                            NDEX_C_EQ (MAX_ACT, N_CEQ),      &
                            I_DIAG (N_D_FRE)
 REAL(DP), INTENT(INOUT) :: S_SKY (N_COEFF), CC (N_D_FRE)
 REAL(DP), INTENT(INOUT) :: S_SKY_LOWER (NOT_SYM)
 INTEGER  :: IC, IEQ, NEQ, NTEST
                                      
! N_D_FRE         = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM
! N_COEFF         = NUMBER OF NON-ZERO COEFFICIENTS IN S_SKY
! NOT_SYM         = NUMBER OF NON-ZERO COEFF IN S_SKY_LOWER, 
!                   = 0 OR N_COEFF
! S_SKY           = SYS. EQ. SQ. MATRIX
! CC              = SYS. EQ. COL. MATRIX
! N_G_DOF         = NO. PARAMETERS PER NODE
! N_RES_EQ        = NO OF CONSTRAINT EQS OF EACH TYPE
! COEF_C_EQ (I,J) = CONSTRAINT EQS COEFF I FOR EQ J
! NDEX_C_EQ (I,J) = CONSTRAINT EQS DOF NO I FOR EQ J
! N_CEQ           = NUMBER OF CONSTRAINT EQUATIONS

  IF ( I_BUG == 1 ) THEN
    PRINT *, 'Entering SKY_BC_GEN'
  END IF
!  WORK BACKWARD FROM LARGEST CONSTRAINT TYPE (DO TYPE 1 LAST)
  IF ( SUM(N_RES_EQ) < 1) THEN
    PRINT *, 'WARNING: NO RESTRAINTS APPLIED'
    N_WARN = N_WARN + 1 ! INCREMENT WARNING
  END IF
  IEQ = SUM ( N_RES_EQ ) + 1  
!b XXX if MAX_ACT > 1, split into 2 passes: LCE1, LCE2
  DO IC = MAX_ACT, 1, -1
    NTEST = N_RES_EQ (IC)
    IF ( NTEST > 0 ) THEN ! Apply constraint of type IC
      IF ( IC == 1 ) THEN
!-->         TYPE 1     D(L1) = C1
        DO NEQ = 1, NTEST
          IEQ = IEQ - 1 
          CALL SKY_TYPE_1_GEN (NDEX_C_EQ (1, IEQ), &
               COEF_C_EQ (1, IEQ), S_SKY, CC, I_DIAG, S_SKY_LOWER)
        END DO
      ELSE
        IF ( .NOT. SYMMETRIC )                         &
           STOP 'NOT SYMMETRIC, SKY_BC_GEN for TYPE > 1'
!-->         TYPE 2   D(L1) + C1*D(L2) = C2
!-->         TYPE 3   D(L1) + C1*D(L2) + C2*D(L3) = C3, ETC.
        DO NEQ = 1, NTEST
          IEQ = IEQ - 1 
          !b CALL SKY_LCE (S_SKY, CC, IC, NDEX_C_EQ (1:IC, IEQ), &
                        !b COEF_C_EQ (1:IC,IEQ), I_DIAG)
          CALL SKY_LCE_SYM (S_SKY, CC, IC, NDEX_C_EQ (1:IC, IEQ), &
                            COEF_C_EQ (1:IC,IEQ), I_DIAG)
        END DO
      END IF
    END IF
  END DO
END SUBROUTINE SKY_BC_GEN
                         
SUBROUTINE CHECK_SKY_SYS (S, C, I_DIAG, S_LOWER) !jea
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! WARNING, this code is not fully checked. No known bugs.
!  CHECK SKYLINE SYSTEM FOR INVALID EQUATIONS & WARN
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 INTEGER,  INTENT(IN)    :: I_DIAG (N_D_FRE)
 REAL(DP), INTENT(INOUT) :: S (N_COEFF), C (N_D_FRE)
 REAL(DP), INTENT(INOUT) :: S_LOWER (NOT_SYM) !jea
 INTEGER  :: I, J, K
 REAL(DP) :: SMAX, TEST, ZERO = 0.d0
 LOGICAL  :: WARN

! C       = SYSTEM COLUMN MATRIX
! I_DIAG  = POINTER TO DIAGONAL COEFFICIENT IN S
! S       = SYSTEM SQUARE MATRIX IN SKY MODE (upper)
! S_LOWER = SYSTEM SQUARE MATRIX IN SKY MODE (lower) !jea
! N_COEFF = NUMBER OF COEFFICIENTS IN S
! NOT_SYM = NUMBER OF COEFFICIENTS IN S_LOWER !jea
! N_D_FRE = NUMBER OF EQUATIONS

  SMAX = ZERO; WARN = .FALSE.    ! INITALIZE
  DO I = 1, N_D_FRE
    !b TEST = ABS (S (I_DIAG (I) ) )
    TEST = S (I_DIAG (I) ) 
    IF ( ABS (TEST) > SMAX) SMAX = TEST
    IF ( TEST <= 0.d0 ) WARN = .TRUE.
  END DO

  IF ( WARN ) THEN
    WRITE (N_PRT, "(/, 'WARNINGS FROM CHECK_SKY_SYS:')")
    N_WARN = N_WARN + 1 ! INCREMENT WARNING
  END IF
  IF ( SMAX <= ZERO ) STOP &
      'ALL ELEMENT STIFFNESSES ZERO, CHECK_SKY_SYS'
  K = 0
  DO I = 1, MAX_NP
    DO J = 1, N_G_DOF
      K = K + 1
      TEST = S (I_DIAG (K) )
      IF ( TEST <= ZERO ) THEN
        IF ( TEST == ZERO ) THEN
          WRITE (N_PRT, 200) I, J
          200 FORMAT ('WARNING, NODE ',I5,' DOF',I3,' RESTRAINED')
          N_WARN = N_WARN + 1 ! INCREMENT WARNING
        END IF
        IF ( TEST < ZERO ) THEN
          WRITE (N_PRT, 300) I, J
          300 FORMAT ('ERROR, NODE ',I5,' DOF',I3,' RESTRAINED')
          N_WARN = N_WARN + 1 ! INCREMENT WARNING
        END IF
!        SET DOF K TO ZERO IN SKYLINE MODE
        CALL SKY_TYPE_1_GEN (K, ZERO, S, C, I_DIAG, S_LOWER) !jea
      END IF
    END DO
  END DO
END SUBROUTINE CHECK_SKY_SYS
                            
SUBROUTINE SKY_DIAG (I_DOF_HI, I_DIAG)
! * * * * * * * * * * * * * * * * * * * * * * * * *
! USE COLUMN HEIGHTS TO FIND DIAGONAL LOCATIONS    
! FOR SKYLINE STORAGE MODE
! * * * * * * * * * * * * * * * * * * * * * * * * *
!  ASSUMING SYMMETRIC COLUMNS STORED FROM TOP DOWN
Use System_Constants
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: I_DOF_HI (N_D_FRE)
 INTEGER, INTENT(OUT) :: I_DIAG (N_D_FRE)
 INTEGER  :: I, IPOINT

! N_D_FRE     = TOTAL NO OF SYSTEM EQUATIONS
! I_DOF_HI(I) = COLUMN HEIGHT OF EQ I, WITH DIAG
! I_DIAG(I)   = LOCATION OF DIAG OF I-TH EQ IN UPPER TRIANGLE

!    TOTAL NUMBER OF SQ MATRIX TERMS = I_DIAG(N_D_FRE)
  IPOINT = 0
  DO I = 1, N_D_FRE
    IPOINT     = IPOINT + I_DOF_HI (I)
    I_DIAG (I) = IPOINT
  END DO
END SUBROUTINE SKY_DIAG
                      
FUNCTION GET_SKY_DIAG (I_DOF_HI) RESULT (I_DIAG)
! * * * * * * * * * * * * * * * * * * * * * * * * *
! USE COLUMN HEIGHTS TO FIND DIAGONAL LOCATIONS    
!           FOR SKYLINE STORAGE MODE
! * * * * * * * * * * * * * * * * * * * * * * * * *
!  ASSUMING SYMMETRIC COLUMNS STORED FROM TOP DOWN
Use System_Constants
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: I_DOF_HI (N_D_FRE)
 INTEGER              :: I_DIAG (N_D_FRE)
 INTEGER  :: I, IPOINT

! N_D_FRE     = TOTAL NO OF SYSTEM EQUATIONS
! I_DOF_HI(I) = COLUMN HEIGHT OF EQ I, WITH DIAG
! I_DIAG(I)   = LOCATION OF DIAG OF I-TH EQ IN UPPER TRIANGLE

!    TOTAL NUMBER OF SQ MATRIX TERMS = I_DIAG(N_D_FRE)
  IPOINT = 0
  DO I = 1, N_D_FRE
    IPOINT     = IPOINT + I_DOF_HI (I)
    I_DIAG (I) = IPOINT
  END DO
END FUNCTION GET_SKY_DIAG
                      
SUBROUTINE SKY_HI (NODES, I_DOF_HI)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! FIND COLUMN HEIGHTS OF SYSTEM EQUATIONS IN SKYLINE 
!    STORAGE MODE (ZERO HEIGHT IS INACTIVE EQ)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data       ! for ELEM_NODES (LT_N), INDEX (LT_FREE)
Use Interface_Header
Use Sys_Properties_Data
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: NODES (L_S_TOT, NOD_PER_EL)
 INTEGER, INTENT(OUT) :: I_DOF_HI (N_D_FRE)
!                  Automatic Arrays
 INTEGER :: LT_NODES (NOD_PER_EL), LT_INDEX (N_EL_FRE) 
 INTEGER :: L_HIGH (N_EL_FRE)
 INTEGER :: IE, J, NDX

! DELETED     = LOGICAL ELEM FLAG, TRUE IF DELETED
! ELEM_NODES  = ELEMENT NODAL INCIDENCES
! INDEX(I)    = SYS DOF NUMBER OF ELEMENT DOF I
! I_DOF_HI(I) = COL HEIGHT OF SYS DOF I
! N_D_FRE     = TOTAL NO OF SYSTEM DOF
! N_EL_FRE    = MAX NUMBER OF ELEMENT PARAMETERS (DOF)
! N_ELEMS     = NUMBER OF ELEMENTS
! NOD_PER_EL  = NUMBER OF NODES PER ELEMENT
! N_G_DOF     = NUMBER OF PARAMETERS PER NODE
! NODES       = NODAL INCIDENCES OF ALL ELEMENTS
 
!  ZERO HEIGHTS
  I_DOF_HI = 0 

!  LOOP OVER ELEMENTS
  DO IE = 1, N_ELEMS ! (boundary segments use same nodes)
    IF ( DELETED (IE) ) CYCLE ! to next element
      LT_NODES = 0; LT_INDEX = 0 !b
!    EXTRACT NODES,  FIND DOF NOS
      LT_N = LT_DATA (1, L_TYPE(IE)) !b CALL ...
      !b LT_NODES = GET_ELEM_NODES (IE, LT_N, NODES)    
      LT_NODES (1:LT_N) = GET_ELEM_NODES (IE, LT_N, NODES)    

!    FIND ELEMENT COLUMN HEIGHTS
      LT_FREE = LT_N*N_G_DOF
      !b LT_INDEX = GET_ELEM_INDEX (LT_N, LT_NODES) 
      LT_INDEX (1:LT_FREE) = GET_ELEM_INDEX (LT_N, LT_NODES) 
      CALL EL_HIGH (LT_INDEX, L_HIGH)
!    COMPARE WITH CURRENT MAXIMUMS
      DO J = 1, LT_FREE
        NDX = LT_INDEX (J)
        IF ( NDX > 0 ) THEN
          IF ( I_DOF_HI (NDX) < L_HIGH (J) ) &
               I_DOF_HI (NDX) = L_HIGH (J)
        END IF
      END DO ! for element dof
  END DO ! for all elements
END SUBROUTINE SKY_HI
                     
FUNCTION GET_SKY_HEIGHTS (NODES) RESULT (I_DOF_HI)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! FIND COLUMN HEIGHTS OF SYSTEM EQUATIONS IN SKYLINE 
!    STORAGE MODE (ZERO HEIGHT IS INACTIVE EQ)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data       ! for ELEM_NODES (LT_N), INDEX (LT_FREE)
Use Interface_Header
Use Sys_Properties_Data
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: NODES (L_S_TOT, NOD_PER_EL)
 INTEGER              :: I_DOF_HI (N_D_FRE)
!                  Automatic Arrays
 INTEGER :: LT_NODES (NOD_PER_EL), LT_INDEX (N_EL_FRE)
 INTEGER :: L_HIGH (N_EL_FRE)
 INTEGER :: IE, J, NDX

! DELETED     = LOGICAL ELEM FLAG, TRUE IF DELETED
!   LT_NODES  = ELEMENT NODAL INCIDENCES
! LT_INDEX(I)    = SYS DOF NUMBER OF ELEMENT DOF I
! I_DOF_HI(I) = COL HEIGHT OF SYS DOF I
! N_D_FRE     = TOTAL NO OF SYSTEM DOF
! N_EL_FRE    = MAX NUMBER OF ELEMENT PARAMETERS (DOF)
! N_ELEMS     = NUMBER OF ELEMENTS
! NOD_PER_EL  = NUMBER OF NODES PER ELEMENT
! N_G_DOF     = NUMBER OF PARAMETERS PER NODE
! NODES       = NODAL INCIDENCES OF ALL ELEMENTS
 
!  ZERO HEIGHTS
  I_DOF_HI = 0 
!  LOOP OVER ELEMENTS
  DO IE = 1, N_ELEMS ! (boundary segments use same nodes)
    IF ( DELETED (IE) ) CYCLE ! to next element
!    EXTRACT NODES,  FIND DOF NOS
      LT_N = LT_DATA (1, L_TYPE(IE))
        LT_NODES = GET_ELEM_NODES (IE, LT_N, NODES) 
!    FIND ELEMENT COLUMN HEIGHTS
      LT_FREE = LT_N*N_G_DOF
      LT_INDEX = GET_ELEM_INDEX (LT_N,   LT_NODES) 
      CALL EL_HIGH (LT_INDEX, L_HIGH)
!    COMPARE WITH CURRENT MAXIMUMS
      DO J = 1, LT_FREE
        NDX = LT_INDEX (J)
        IF ( NDX > 0 ) THEN
          IF ( I_DOF_HI (NDX) < L_HIGH (J) ) &
               I_DOF_HI (NDX) = L_HIGH (J)
        END IF
      END DO ! for element dof
  END DO ! for all elements
END FUNCTION GET_SKY_HEIGHTS
                     
SUBROUTINE SKY_LCE (S, C, NCD, NDX, A, I_DIAG)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! APPLY LINEAR CONSTRAINTS TO SKYLINE SYMMETRIC EQUATIONS
!             S*D = C,  WITH S IN VECTOR MODE AND
! D(NDX(1))+A(1)*D(NDX(2))+...A(NCD-1)*D(NDX(NCD))=A(NCD)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 REAL(DP), PARAMETER     :: ZERO = 0.d0, ONE = 1.d0
 INTEGER,  INTENT(IN)    :: NCD, I_DIAG (N_D_FRE), NDX (NCD)
 REAL(DP), INTENT(IN)    :: A (NCD)
 REAL(DP), INTENT(INOUT) :: S (N_COEFF), C (N_D_FRE)

 LOGICAL  :: SKIP
 INTEGER  :: I, IBOT, I_J_V, I_NR_V, ITOP, J, K, L, LL, NR
 REAL(DP) :: CR, E, ESC, SRR, SRRP1
!
! N_COEFF = NUMBER OF TERMS FOR SQ MATRIX IN VECTOR MODE
! NCD     = TOTAL NUMBER OF DOF IN CONSTRAINT EQUATION
! N_D_FRE = TOTAL NUMBER OF DEGREES OF FREEDOM
! NDX(I)  = SYS DOF NOS OF CONSTRAINT TERM I
! NR      = REDUNDANT DEGREE OF FREEDOM = NDX(1)
! A(J)    = NORMALIZED COEFF OF (J+1) TERM, A0 = 1.0

!                  S*D=C SYSTEM PARTITIONED
! C = SYSTEM COLUMN MATRIX  | SII | SIR | SID | |DI| |CI|
! S = SYSTEM SQUARE MATRIX  |-----|-----|-----| |--| |--|
! I-INDEPENDENT             | SRI | SRR | SRD |*|DR|=|CR|
! R-REDUNDANT               |-----|-----|-----| |--| |--|
! D-DEPENDENT               | SDI | SDR | SDD | |DD| |CD|
 
  IF ( I_BUG > 0 ) PRINT *, 'Entering SKY_LCE'
  IF ( ANY (NDX < 1)       ) STOP 'SKY_LCE: BAD CONSTRAINT INDEX'
  IF ( ANY (NDX > N_D_FRE) ) STOP 'SKY_LCE: BAD CONSTRAINT INDEX'
  E     = A (NCD)
  NR    = NDX (1)
  SRR   = S (I_DIAG (NR) )
  SRRP1 = SRR + ONE
  CR    = C (NR)
  ESC   = E * SRRP1 - CR
!  FIND TOP AND BOTTOM OF THE SKYLINE OF THE NR COLUMN
  CALL SKY_TOP_BOT (I_DIAG, NR, ITOP, IBOT)
!  FORM MODIFIED COLUMN MATRIX, CX = CX - E*SXR
  IF ( E /= ZERO ) THEN
    DO I = ITOP, IBOT
      !b IF ( I_BUG > 0 ) PRINT *, 'I, NR ', I, NR
      CALL SUB_IN_SKY (I_DIAG, I, NR, I_NR_V)
      IF ( I_NR_V /= 0 ) C (I) = C (I) - E * S (I_NR_V)
    END DO
  END IF
!  ADDITIONAL COLUMN CHANGES FOR CD AND CR
  IF ( NCD > 1 .AND. ESC /= ZERO ) THEN
    DO K = 2, NCD
      C (NDX (K) ) = C (NDX (K) ) + A (K - 1) * ESC
    END DO
  END IF
  C (NR) = E
!    *** SQUARE MATRIX COLUMN MODIFICATIONS ***
  IF ( NCD > 1)  THEN
!    FORM SID,  BEGIN SDD
    DO K = 2, NCD
      J = NDX (K)
      DO I = ITOP, IBOT
        !b IF ( I_BUG > 0 ) PRINT *, 'I, J  ', I, J 
        !b IF ( I_BUG > 0 ) PRINT *, 'I, NR ', I, NR
        CALL SUB_IN_SKY (I_DIAG, I, J, I_J_V)
        CALL SUB_IN_SKY (I_DIAG, I, NR, I_NR_V)
        IF ( I_NR_V /= 0 ) THEN 
          IF ( S (I_NR_V) /= ZERO ) THEN 
            IF ( I_J_V == 0 ) STOP 'INVALID ADDRESS, SKY_LCE'
!         DO  NOT MODIFY NR  COLUMN, MAY NEED IT LATER
            SKIP = .FALSE.
            DO LL = K + 1, NCD
              IF ( I == NDX (LL) ) SKIP = .TRUE.
            END DO
            IF ( (.NOT. SKIP) .AND. I /= NR ) &
              S (I_J_V) = S (I_J_V) - S (I_NR_V) * A (K - 1)
          END IF ! S (I_NR_V) 
        END IF ! I_NR_V 
      END DO
    END DO
!  COMPLETE SDD
    DO K = 2, NCD
      I = NDX (K)
      DO L = 2, K
        J = NDX (L)
        !b IF ( I_BUG > 0 ) PRINT *, 'I, J  ', I, J 
        !b IF ( I_BUG > 0 ) PRINT *, 'I, NR ', I, NR
        CALL SUB_IN_SKY (I_DIAG, I, J, I_J_V)
        CALL SUB_IN_SKY (I_DIAG, I, NR, I_NR_V)
        IF ( I_J_V == 0 .OR. I_NR_V == 0 )  &
             STOP 'INVALID ADDRESS, SKY_LCE'
        S (I_J_V) = S (I_J_V) + SRRP1 * A (K - 1) * A (L - 1) &
                  - A (L - 1) * S (I_NR_V)
      END DO
    END DO
  END IF
!b XXX end LCE1 begin LCE2, split into 2 routines
!        *** INSERT CONSTRAINT EQUATION ***
! WARNING: NEXT LOOP NOT VALID FOR COUPLED CONSTRAINTS
  DO I = ITOP, IBOT
    !b IF ( I_BUG > 0 ) PRINT *, 'I, NR ', I, NR
    CALL SUB_IN_SKY (I_DIAG, I, NR, I_NR_V)
    IF ( I_NR_V /= 0 ) S (I_NR_V) = ZERO
  END DO
  S (I_DIAG (NR)) = 1.d0
  IF ( NCD > 1 ) THEN
    DO K = 2, NCD
      I = NDX (K)
      !b IF ( I_BUG > 0 ) PRINT *, 'I, NR ', I, NR
      CALL SUB_IN_SKY (I_DIAG, I, NR, I_NR_V)
      IF ( I_NR_V == 0 ) STOP 'INVALID ADDRESS, SKY_LCE'
      S (I_NR_V) = A (K - 1)
    END DO
!   *** MODIFICATIONS COMPLETED, CHECK DIAGONAL ***
    DO K = 2, NCD
      I = NDX (K)
      IF ( S (I_DIAG (I)) <= ZERO ) THEN ; WRITE (6,  * ) &
        'WARNING: NEGATIVE DIAGONAL FOR CONSTRAINT SET', NDX
        N_WARN = N_WARN + 1 ! INCREMENT WARNING
      END IF
    END DO
  END IF
  IF ( I_BUG > 0 ) PRINT *, 'Exiting  SKY_LCE'
END SUBROUTINE SKY_LCE
                      
SUBROUTINE SKY_SOLVE_SYM (A, B, X, I_DIAG, FACT, BACK)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! PERFORM  A =(U)T*D*U FACTORIZATION AND/OR BACKSUBSTITUTION
! OF SYMMETRIC POSITIVE DEFINITE SYSTEM OF EQUATIONS A*X = B
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants     ! for N_COEFF, N_D_FRE
 IMPLICIT NONE
              
 LOGICAL,  INTENT(IN)    :: FACT, BACK
 INTEGER,  INTENT(IN)    :: I_DIAG (N_D_FRE)
 REAL(DP), INTENT(INOUT) :: A (N_COEFF), B (N_D_FRE)
 REAL(DP), INTENT(OUT)   :: X (N_D_FRE)
                                       
 INTEGER  :: I, ID, IE, IH, IR, IS, J, JD, JH, JR, K
 REAL(DP) :: D, ENERGY, ZERO = 0.d0
!
! A(N_COEFF)      = UPPER TRIANGULAR COEFFICIENT MATRIX STORED IN
!                   COLUMN FORM (HOLDS D & U ON EXIT)
! B(N_D_FRE)      = RIGHT SIDE VECTOR (HOLDS X ON EXIT)
! X(N_D_FRE)      = SOLUTION VECTOR
! I_DIAG(N_D_FRE) = ADDRESSES OF DIAGONAL TERMS IN A(N_COEFF)
! N_D_FRE         = NUMBER OF EQUATIONS
! FACT            = .TRUE., FACTOR A(N_COEFF)
!                   .FALSE., DO NOT FACTOR A(N_COEFF)
! BACK            = .TRUE., FORWARD REDUCE B(N_D_FRE) & BACKSUBSTITUTE
!                   .FALSE., DO NOT FORWARD REDUCE & BACKSUBSTITUTE
!
!.... FACTOR A,  REDUCE B
!
  JR = 0
  DO J = 1, N_D_FRE
    JD = I_DIAG (J)
    JH = JD-JR
    IS = J - JH + 2
    IF ( JH <  2 ) GO TO 600
    IF ( JH == 2 ) GO TO 300
100 IF ( .NOT. FACT ) GO TO 500
    IE = J - 1
    K = JR + 2
    ID = I_DIAG (IS - 1)
!
!.... REDUCE ALL EQUATIONS EXCEPT DIAGONAL
!
    DO I = IS, IE
      IR = ID
      ID = I_DIAG (I)
      IH = MIN0 (ID-IR - 1, I - IS + 1)
      IF ( IH > 0 ) A (K) = A (K) &
                    - DOT_PRODUCT (A (K-IH:K-1), A (ID-IH:ID-1))
      K = K + 1
    END DO ! on I
!
!.... REDUCE DIAGONAL TERM
!
300 IF ( .NOT. FACT ) GO TO 500
    IR = JR + 1
    IE = JD-1
    K = J - JD
    DO I = IR, IE
      ID = I_DIAG (K + I)
      IF ( A (ID) == ZERO ) THEN
        STOP 'ZERO DIAGONAL IN SKY_SOLVE_SYM'
      ELSE
        D = A (I)
        A (I) = A (I) / A (ID)
        A (JD) = A (JD) - D * A (I)
      END IF
    END DO ! on I
!
!.... REDUCE RHS
!
500 IF ( BACK ) B (J) = B (J) &
                - DOT_PRODUCT(A (JR+1:JR+JH-1), B (IS-1:IS+JH-3))
600 CONTINUE
    JR = JD
  END DO ! on J
  IF ( .NOT. BACK ) RETURN
!
!.... DIVIDE BY DIAGONAL PIVOTS
!
  ENERGY = 0.d0 !b
  DO I = 1, N_D_FRE
    D  = B (I)
    ID = I_DIAG (I)
    IF ( A (ID) /= ZERO ) B (I) = B (I) / A (ID)
    ENERGY = ENERGY + D * B (I)
  END DO
  !b PRINT *, ' '
  !b PRINT *, 'EQUATION STRAIN ENERGY VALUE = ', ENERGY * 0.5d0 
  IF ( ENERGY < 0.d0 ) THEN
    PRINT *, 'WARNING, SKY_SOLVE_SYM: NEGATIVE ENERGY IN SYSTEM'
    PRINT *, 'WARNING: IS unsymmetric NEEDED AS KEYWORD INPUT ?'
    N_WARN = N_WARN + 2
    PRINT *, 'ERROR, SKY_SOLVE_SYM: SYMMETRIC SOLVER FATAL ERROR'     
    STOP     'SKY_SOLVE_SYM: FATAL ERROR, NEGATIVE ENERGY'
  ELSE
    !b PRINT *, 'EQUATION STRAIN ENERGY NORM  = ', SQRT (ENERGY)
  END IF ! INVALID SYSTEM
!
!.... BACKSUBSTITUTE
!
  J = N_D_FRE
  JD = I_DIAG (J)
800 D = B (J)  !9 replace with while on J
  J = J - 1
  IF ( J <= 0 ) RETURN
  JR = I_DIAG (J)
  IF ( JD-JR > 1 ) THEN 
    IS = J - JD+JR + 2
    K = JR - IS + 1
    DO I = IS, J
      B (I) = B (I) - A (I + K) * D
    END DO
    X = B
  END IF
  JD = JR
  GO TO 800
END SUBROUTINE SKY_SOLVE_SYM
                            
SUBROUTINE SKY_FORWARD_SYM (A, B, I_DIAG)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! PREFORM FORWARD SUBSTITUTION, ONLY, FOR PREVIOUSLY
! FACTORED SYMMETRIC SKYLINE MATRIX SYSTEM, A*X = B VIA
! SKY_SOLVE_SYM (A, B, X, I_DIAG, .TRUE., .FALSE.)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants     ! for N_COEFF, N_D_FRE
 IMPLICIT NONE

 INTEGER,  INTENT(IN)    :: I_DIAG (N_D_FRE)
 REAL(DP), INTENT(IN)    :: A (N_COEFF)
 REAL(DP), INTENT(INOUT) :: B (N_D_FRE)
 !b REAL(DP), INTENT(IN)    :: X (N_D_FRE)

 INTEGER  :: IS, J, JD, JH, JR

! A(N_COEFF)      = UPPER TRIANGULAR COEFFICIENT MATRIX STORED IN
!                   COLUMN FORM (HOLDS D & U ON EXIT)
! B(N_D_FRE)      = RIGHT SIDE VECTOR (HOLDS X ON EXIT)
! X(N_D_FRE)      = SOLUTION VECTOR
! I_DIAG(N_D_FRE) = ADDRESSES OF DIAGONAL TERMS IN A(N_COEFF)
! N_D_FRE         = NUMBER OF EQUATIONS

  JR = 0
  DO J = 1, N_D_FRE
    JD = I_DIAG (J)
    JH = JD-JR
    IS = J - JH + 2
    !b IF ( JH <  2 ) GO TO 600
    IF ( JH > 1 ) THEN      
      B (J) = B (J) - DOT_PRODUCT(A (JR+1:JR+JH-1), B (IS-1:IS+JH-3))
    END IF
!b 600 CONTINUE
    JR = JD
  END DO ! on J
END SUBROUTINE SKY_FORWARD_SYM 

SUBROUTINE SKY_SOLVE_UNSYM (A, C, B, X, I_DIAG, FACT, BACK)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! PROGRAM TO PERFORM (C\A)=L*D*U FACTORIZATION AND/OR BACKSUBSTITUTION
!    OF AN UNSYMMETRIC SYSTEM OF EQUATIONS, (C\A)*X = B
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants     ! for N_COEFF, N_D_FRE
 IMPLICIT NONE
              
 LOGICAL,  INTENT(IN)    :: FACT, BACK
 INTEGER,  INTENT(IN)    :: I_DIAG (N_D_FRE)
 REAL(DP), INTENT(INOUT) :: A (N_COEFF), C (N_COEFF), B (N_D_FRE)
 REAL(DP), INTENT(OUT)   :: X (N_D_FRE)
                                       
 INTEGER  :: I, ID, IR, J, JD, JH, JR, K, M, N_END, N_START, NT
 REAL(DP) :: D
!
! A(N_COEFF)      = UPPER TRIANGULAR COEFFICIENT MATRIX STORED IN
!                   COLUMN FORM (HOLDS D & U ON EXIT)
! C(N_COEFF)      = LOWER TRIANGULAR COEFFICIENT MATRIX STORED IN
!                   COLUMN FORM (HOLDS L ON EXIT), DIAGONAL IGNORED
! B(N_D_FRE)      = RIGHT SIDE VECTOR (COULD HOLD X ON EXIT)
! X(N_D_FRE)      = SOLUTION VECTOR
! I_DIAG(N_D_FRE) = ADDRESSES OF DIAGONAL TERMS IN A(N_COEFF) AND C
! N_COEFF         = NUMBER OF STORED TERMS IN A AND C, = I_DIAG(N_D_FRE)
! N_D_FRE         = NUMBER OF EQUATIONS
! FACT            = .TRUE., FACTOR A(N_COEFF) AND C(N_COEFF)
!                   .FALSE., DO NOT FACTOR A(N_COEFF) AND C(N_COEFF)
! BACK            = .TRUE., FORWARD REDUCE B(N_D_FRE) & BACKSUBSTITUTE
!                   .FALSE., DO NOT FORWARD REDUCE & BACKSUBSTITUTE
!
  JR = 0
  DO 15 J = 1, N_D_FRE
    JD = I_DIAG (J)
    JH = JD-JR
    IF ( JH <= 1 ) GO TO 10
    N_START = J + 1 - JH
    N_END = J - 1
    IF ( .NOT. FACT ) GO TO 20
    K = JR + 1
    ID = 0
!
!.... REDUCE ALL EQUATIONS EXCEPT DIAGONAL
!
    DO 30 I = N_START, N_END
      IR = ID
      ID = I_DIAG (I)
      NT = MIN0 (ID-IR - 1, I - N_START)
      IF ( NT == 0 ) GO TO 40
      A (K) = A (K) - DOT_PRODUCT (A (K-NT:K-1), C (ID-NT:ID-1))
      C (K) = C (K) - DOT_PRODUCT (C (K-NT:K-1), A (ID-NT:ID-1))
40    IF ( A (ID) /= 0. ) THEN
        C (K) = C (K) / A (ID)
      ELSE
        STOP 'ZERO DIAGONAL IN SKY_SOLVE_UNSYM'
      END IF
      K = K + 1
30  END DO
!
!.... REDUCE DIAGONAL TERM
!
    A (JD) = A (JD) &
           - DOT_PRODUCT (A (JR+1:JR+JH-1), C (JR+1:JR+JH-1))
!
!.... FORWARD REDUCE THE R.H.S.
!
20  IF ( BACK ) B (J) = B (J) &
      - DOT_PRODUCT (C (JR+1:JR+JH-1), B (N_START:N_START+JH-2))
10  JR = JD
15 END DO
  IF ( .NOT. BACK ) RETURN
!
!.... BACKSUBSTITUTION
!
  J = N_D_FRE
  JD = I_DIAG (J)
  IF ( A (JD) /= 0. ) B (J) = B (J) / A (JD)
  IF ( N_D_FRE == 1 ) THEN
     X = B ; RETURN
  END IF
50 D = B (J)
  J = J - 1
  JR = I_DIAG (J)
  IF ( JD-JR <= 1 ) GO TO 60
  M = J - JD+JR + 2
  K = JR - M + 1
  DO 70 I = M, J
    B (I) = B (I) - A (I + K) * D
70 END DO
60 JD = JR
  IF ( A (JD) /= 0. ) B (J) = B (J) / A (JD)
  IF ( J > 1 ) GO TO 50
!   MOVE INTO X
  X = B
END SUBROUTINE SKY_SOLVE_UNSYM
                              
SUBROUTINE SKY_STORE_SQ (LT_FREE, INDEX, I_DIAG, S, SS, SS_LOWER)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!   STORE ELEMENT SQUARE MATRIX TO SYSTEM SQUARE MATRIX STORED
!   IN SYMMETRIC OR NON-SYMMETRIC SKYLINE MODE (Revised 5/9/98)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 INTEGER,  INTENT(IN)    :: LT_FREE
 INTEGER,  INTENT(IN)    :: INDEX (LT_FREE), I_DIAG (N_D_FRE)
 REAL(DP), INTENT(INOUT) :: SS (N_COEFF), SS_LOWER (NOT_SYM)
 REAL(DP), INTENT(IN)    :: S (LT_FREE, LT_FREE)
 INTEGER  :: I, J, NDX_I, NDX_J, NDX_V
                                                    
! N_COEFF   = NO COEFF IN UPPER HALF OF SQ MATRIX = I_DIAG(N_D_FRE)
! NOT_SYM   = NO COEFF IN LOWER HALF OF SQ MATRIX = I_DIAG(N_D_FRE)
! N_D_FRE   = TOTAL NO OF DOF IN SYSTEM
! LT_FREE   = NUMBER OF ELEMENT DEGREES OF FREEDOM
! INDEX(I)  = SYS DOF NO OF ELEMENT DOF I
! I_DIAG(I) = LOCATION OF DIAGONAL OF I-TH EQ
! S         = ELEMENT SQUARE MATRIX
! SS        = UPPER HALF OF SYS SQ MATRIX IN SKYLINE VECTOR MODE
! SS_LOWER  = LOWER HALF OF SYS SQ MATRIX IN SKYLINE VECTOR MODE
! SYMMETRIC = .true. if symmetric so: ss_lower = 0, NOT_SYM = 0
                                                               
!  LOOP OVER ELEMENT COEFFICIENTS
  DO J = 1, LT_FREE                    ! element columns
    NDX_J = INDEX (J)                  ! system column
!    ALLOW FOR OMITTED NODES
    IF ( NDX_J > 0 ) THEN              ! equation exists
      DO I = 1, LT_FREE                ! element rows
        NDX_I = INDEX (I)              ! system row
        IF ( NDX_I > 0 ) THEN          ! equation exists
          IF ( I_BUG > 0 .AND. SYMMETRIC .AND. I == J &
            .AND. S (I, J) < 0.d0) THEN
            PRINT *,'BAD DIAGONAL ELEMENT TERM = ', S (I, J)
            PRINT *,'ELEMENT ROW = ', I
            PRINT *,'SYSTEM  ROW = ', NDX_I
            STOP    'NEGATIVE ELEMENT DIAGONAL, SKY_STORE_SQ'
          END IF
!          FIND SYSTEM COEFF IN VECTOR SS (AND SS_LOWER)
          NDX_V = I_DIAG ( MAX0 (NDX_I, NDX_J) ) - IABS (NDX_J-NDX_I)
          IF ( NDX_I <= NDX_J ) THEN   ! Upper half
            SS (NDX_V) = SS (NDX_V) + S (I, J)
          ELSE                         ! Lower half
            IF ( .NOT. SYMMETRIC ) &
              SS_LOWER (NDX_V) = SS_LOWER (NDX_V) + S (I, J)
          END IF
        END IF
      END DO
    END IF
  END DO
END SUBROUTINE SKY_STORE_SQ
                           
SUBROUTINE SKY_TOP_BOT (I_DIAG, J_COL, J_TOP, J_BOT)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! FIND THE TOP (J_TOP) AS WELL AS (MAXIMUM) BOTTOM
! (J_BOT) INDEX  OF THE SKYLINE OF COL "J_COL"
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! ASSUMING SYMM EQS, COLS STORED FROM TOP DOWN
Use System_Constants
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: I_DIAG (N_D_FRE), J_COL
 INTEGER, INTENT(OUT) :: J_TOP, J_BOT
 INTEGER  :: I, MINI
 LOGICAL  :: DEBUG = .FALSE.

! N_D_FRE   = TOTAL NO OF SYSTEM EQUATIONS
! I_DIAG(I) = LOCATION OF DIAG OF I-TH EQ
!
! * FINDING J_TOP
  IF ( J_COL > 1 ) THEN
    J_TOP = J_COL - (I_DIAG (J_COL) - I_DIAG (J_COL - 1) ) + 1
  ELSE
    J_TOP = 1
  END IF
        
! * FINDING J_BOT
  J_BOT = J_COL
  DO I = J_COL + 1, N_D_FRE
    MINI = I - (I_DIAG (I) - I_DIAG (I - 1) ) + 1
    IF ( MINI <= J_COL ) J_BOT = I
  END DO
  IF ( DEBUG ) WRITE (N_BUG, *) "SKY_TOP_BOT: ", J_COL, J_TOP, J_BOT
END SUBROUTINE SKY_TOP_BOT
                          
SUBROUTINE SUB_IN_SKY (I_DIAG, I, J, I_J_V)
! * * * * * * * * * * * * * * * * * * * * * * *
! CONVERT (I,J) FULL SYMMETRIC MATRIX SUBSCRIPTS TO
! (I_J_V) SUBSCRIPT OF VECTOR SKYLINE STORAGE MODE.
! I_J_V = 0 IF OUTSIDE SKYLINE AND THUS HAS 0 VALUE.
! * * * * * * * * * * * * * * * * * * * * * * *
! ASSUMING SYMM EQS, COLS STORED FROM TOP DOWN
! REVISED FOR I_J_V = 0, 9/97
Use System_Constants
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: I_DIAG (N_D_FRE), I, J
 INTEGER, INTENT(OUT) :: I_J_V
 INTEGER  :: ID

! N_D_FRE   = TOTAL NO OF SYSTEM EQUATIONS
! I_DIAG(I) = LOCATION OF DIAG OF I-TH EQ
  ID  = MAX0 (I,J)
  !b if ( i < 1 .or. I > N_D_FRE ) stop 'invalid diag. SUB_IN_SKY'
  !b if ( j < 1 .or. j > N_D_FRE ) stop 'invalid diag. SUB_IN_SKY'
  I_J_V = I_DIAG (ID) - IABS(I-J)
! test for out of bounds
  if (I_J_V < 1 .or. I_J_V > n_coeff ) then
    print *, 'out of bounds, SUB_IN_SKY: I, J, V = ', I, J, I_J_V
    CALL SKY_FUL_LOCATIONS (I_DIAG)          ! locations
    stop 'out of bounds, SUB_IN_SKY'
  end if
!   TEST FOR ZERO ENTRY OUTSIDE SKYLINE
  IF ( ID > 1 ) THEN
    IF ( I_J_V <= I_DIAG (ID - 1) ) I_J_V = 0
  END IF
END SUBROUTINE SUB_IN_SKY
                         
FUNCTION GET_SKY_SUBSCRIPT (I_DIAG, I, J) RESULT (I_J_V)
! * * * * * * * * * * * * * * * * * * * * * * *
! CONVERT (I,J) FULL SYMMETRIC MATRIX SUBSCRIPTS TO
! (I_J_V) SUBSCRIPT OF VECTOR SKYLINE STORAGE MODE.
! I_J_V = 0 IF OUTSIDE SKYLINE AND THUS HAS 0 VALUE.
! * * * * * * * * * * * * * * * * * * * * * * *
! ASSUMING SYMM EQS, COLS STORED FROM TOP DOWN
! REVISED FOR I_J_V = 0, 9/97
Use System_Constants
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: I_DIAG (N_D_FRE), I, J
 INTEGER              :: I_J_V
 INTEGER  :: ID

! N_D_FRE   = TOTAL NO OF SYSTEM EQUATIONS
! I_DIAG(I) = LOCATION OF DIAG OF I-TH EQ
  ID  = MAX0 (I,J)
  IF ( I < 1 .OR. I > N_D_FRE ) STOP 'INVALID INPUT GET_SKY_SUBSCRIPT'
  IF ( J < 1 .OR. J > N_D_FRE ) STOP 'INVALID INPUT GET_SKY_SUBSCRIPT'
  I_J_V = I_DIAG (ID) - IABS(I-J)
! test for out of bounds
  if (I_J_V < 1 .or. I_J_V > n_coeff ) then
    print *, 'out of bounds, GET_SKY_SUBSCRIPT: I, J, V = ', I, J, I_J_V
    CALL SKY_FUL_LOCATIONS (I_DIAG)          ! locations
    stop 'out of bounds, GET_SKY_SUBSCRIPT'
  end if
!   TEST FOR ZERO ENTRY OUTSIDE SKYLINE
  IF ( ID > 1 ) THEN
    IF ( I_J_V <= I_DIAG (ID - 1) ) I_J_V = 0
  END IF
END FUNCTION GET_SKY_SUBSCRIPT
                         
SUBROUTINE CHECK_SKY_SYS_SYM (S, C, I_DIAG)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! WARNING, this code is not fully checked. No known bugs.
!  CHECK SKYLINE SYSTEM FOR INVALID EQUATIONS & WARN
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 REAL(DP), INTENT(INOUT) :: S (N_COEFF), C (N_D_FRE)
 INTEGER,  INTENT(IN)    :: I_DIAG (N_D_FRE)
 INTEGER  :: I, J, K
 REAL(DP) :: SMAX, TEST, ZERO = 0.d0

! C       = SYSTEM COLUMN MATRIX
! I_DIAG  = POINTER TO DIAGONAL COEFFICIENT IN S
! S       = SYSTEM SQUARE MATRIX IN SKYLINE MODE
! N_COEFF = NUMBER OF COEFFICIENTS IN S
! N_D_FRE = NUMBER OF EQUATIONS

  SMAX = ZERO
  DO I = 1, N_D_FRE
    TEST = ABS (S (I_DIAG (I) ) )
    IF ( TEST > SMAX) SMAX = TEST
  END DO
  IF ( SMAX <= ZERO ) STOP &
     'ALL ELEMENT STIFFNESSES ZERO, CHECK_SKY_SYS'
  K = 0
  DO I = 1, MAX_NP
    DO J = 1, N_G_DOF
      K = K + 1
      TEST = S (I_DIAG (K) )
      IF ( TEST <= ZERO ) THEN
        IF ( TEST == ZERO ) THEN ; WRITE (N_PRT, 2) I, J
          2 FORMAT ('WARNING, NODE ',I5,' DOF',I3,' WAS RESTRAINED')
          N_WARN = N_WARN + 1 ! INCREMENT WARNING
        END IF
        IF ( TEST < ZERO ) THEN
          WRITE (N_PRT, 3) I, J
          3 FORMAT ('ERROR, NODE ',I5,' DOF',I3,' WAS RESTRAINED')
          STOP 'ERROR, NEGATIVE DIAGONAL. CHECK_SKY_SYS_SYM'
        END IF ! negative diagonal
!        SET DOF K TO ZERO IN SKYLINE MODE
        CALL SKY_TYPE_1_SYM (K, ZERO, S, C, I_DIAG)
      ENDIF
    END DO
  END DO
END SUBROUTINE CHECK_SKY_SYS_SYM
                                
SUBROUTINE SAVE_SKY_ROWS_SYM (I_BC, S_SKY, CC, I_DIAG)
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
!-->  SAVE INDEPENDENT REACTION EQUATIONS FROM SKYLINE MATRIX
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Interface_Header
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: S_SKY (N_COEFF), CC (N_D_FRE)
 INTEGER,  INTENT(IN) :: I_BC (MAX_NP), I_DIAG (N_D_FRE)
!          Automatic Arrays
 INTEGER :: INDEX (N_G_DOF), KODES (N_G_DOF)
 INTEGER :: J, JKV, J1, J2, IG, INDX, K

! N_G_DOF = NUMBER OF DOF PER NODE
! N_REACT = SEQUENTIAL UNIT TO STORE REACTION DATA
! N_D_FRE = TOTAL NUMBER OF EQUATIONS

  IF ( N_REACT > 0 ) THEN ; REWIND (N_REACT)
  ELSE ; STOP 'NO UNIT FOR REACTION STORAGE, SAVE_SKY_ROWS_SYM'
  END IF
  IF ( I_BUG > 0 ) PRINT *, 'Entering SAVE_SKY_ROWS_SYM'
!  LOOP OVER EQUATIONS FOR ESSENTIAL BC FLAG
  DO J = 1, MAX_NP
    IF ( I_BC (J)  > 0 ) THEN
      CALL PT_CODE (J, I_BC (J), KODES)
      INDEX = GET_INDEX_AT_PT (J) 
      DO IG = 1, N_G_DOF
        IF ( KODES (IG)  > 0 ) THEN
          INDX = INDEX (IG)
!           COLUMN NUMBERS
          J1 = INDX
          J2 = INDX
          IF (INDX > 1) J1 = J2 - (I_DIAG(INDX) - I_DIAG(INDX-1)) + 1
          IF ( INDX < N_D_FRE ) THEN ! FIND LAST COLUMN
            DO K = N_D_FRE, INDX, -1
              CALL SUB_IN_SKY (I_DIAG, INDX, K, JKV)
              IF ( JKV > 0 ) THEN ! FOUND NON-ZERO ENTRY
                J2 = K ; EXIT ! THIS LOOP
              END IF
            END DO
          END IF
!          WRITE NODE, PARAMETER, RANGE OF NON-ZEROS
if ( j1 < 1 ) then !b
  j1 = 1 !b
  N_WARN = N_WARN + 1
  print *, 'WARNING: fixed bounds error, SAVE_SKY_ROWS_SYM' !b
end if !b
if ( j2 > N_D_FRE ) then !b
  j2 = N_D_FRE !b
  N_WARN = N_WARN + 1
  print *, 'fWARNING: ixed bounds error, SAVE_SKY_ROWS_SYM' !b
end if !b
          WRITE (N_REACT) J, IG, J1, J2
!          SAVE ROW OF EQUILIBRIUM EQ
          DO K = J1, J2
            CALL SUB_IN_SKY (I_DIAG, INDX, K, JKV)
              IF ( JKV > 0 ) THEN
                WRITE (N_REACT) S_SKY (JKV)
              ELSE
                WRITE (N_REACT) 0.d0  
              ENDIF
          END DO
          WRITE (N_REACT) CC (INDX)
        ENDIF
      END DO
    ENDIF
  END DO
END SUBROUTINE SAVE_SKY_ROWS_SYM
                                
SUBROUTINE SKY_BC_SYM (N_RES_EQ, COEF_C_EQ, NDEX_C_EQ, S_SKY, &
                       CC, I_DIAG)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! APPLY BOUNDARY CONSTRAINT EQUATIONS TO SYMMETRIC
! SKYLINE SYSTEM OF EQS, S_SKY*D = CC
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 INTEGER,  INTENT(IN)    :: N_RES_EQ  (MAX_ACT)
 INTEGER,  INTENT(IN)    :: NDEX_C_EQ (MAX_ACT, N_CEQ)
 INTEGER,  INTENT(IN)    :: I_DIAG     (N_D_FRE)
 REAL(DP), INTENT(IN)    :: COEF_C_EQ (MAX_ACT, N_CEQ)
 REAL(DP), INTENT(INOUT) :: S_SKY (N_COEFF), CC (N_D_FRE)
 INTEGER  :: IC, IEQ, NEQ, NTEST

! N_D_FRE         = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM
! N_COEFF         = NUMBER OF NON-ZERO COEFFICIENTS IN S_SKY
! S_SKY           = SYS. EQ. SQ. MATRIX
! CC              = SYS. EQ. COL. MATRIX
! N_G_DOF         = NO. PARAMETERS PER NODE
! N_RES_EQ        = NO OF CONSTRAINT EQS OF EACH TYPE
! COEF_C_EQ (I,J) = CONSTRAINT EQS COEFF I FOR EQ J
! NDEX_C_EQ (I,J) = CONSTRAINT EQS DOF NO I FOR EQ J
! N_CEQ           = NUMBER OF CONSTRAINT EQUATIONS

  IF ( I_BUG > 0 ) PRINT *, 'Entering SKY_BC_SYM'
  IEQ = 0
!b XXX if MAX_ACT > 1, split into 2 passes: LCE1, LCE2
!  DO TYPE ONE LAST
  DO IC = MAX_ACT, 1, -1
    NTEST = N_RES_EQ (IC)
    IF ( NTEST > 0 ) THEN
      IF ( IC == 1 ) THEN
!-->         TYPE 1     D(L1) = C1
        DO NEQ = 1, NTEST
          IEQ = IEQ + 1
          CALL SKY_TYPE_1_SYM (NDEX_C_EQ (1, IEQ),         &
                      COEF_C_EQ (1, IEQ), S_SKY, CC, I_DIAG)
        END DO
      ELSE
!-->         TYPE 2   D(L1) + C1*D(L2)            = C2
!-->         TYPE 3   D(L1) + C1*D(L2) + C2*D(L3) = C3, etc
        DO NEQ = 1, NTEST
          IEQ = IEQ + 1
          CALL SKY_LCE_SYM (S_SKY, CC, IC, NDEX_C_EQ (1:IC, IEQ), &
                            COEF_C_EQ (1:IC, IEQ), I_DIAG)
        END DO
      ENDIF  ! SPECIAL TYPE 1, OR GENERAL TYPE
    ENDIF ! TYPE ANY TYPE IC ARE ACTIVE
  END DO
END SUBROUTINE SKY_BC_SYM
                         
SUBROUTINE SKY_LCE_SYM (S, C, NCD, NDX, A, I_DIAG)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! APPLY LINEAR CONSTRAINTS TO SKYLINE SYMMETRIC EQUATIONS
!             S*D = C,  WITH S IN VECTOR MODE AND
! D(NDX(1))+A(1)*D(NDX(2))+...A(NCD-1)*D(NDX(NCD))=A(NCD)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 REAL(DP), PARAMETER     :: ZERO = 0.d0, ONE = 1.d0
 INTEGER,  INTENT(IN)    :: NCD, I_DIAG (N_D_FRE), NDX (NCD)
 REAL(DP), INTENT(IN)    :: A (NCD)
 REAL(DP), INTENT(INOUT) :: S (N_COEFF), C (N_D_FRE)
 LOGICAL  :: SKIP
 INTEGER  :: I, IBOT, I_J_V, I_NR_V, ITOP, J, K, L, LL, NR
 REAL(DP) :: CR, E, ESC, SRR, SRRP1

! N_COEFF = NUMBER OF TERMS FOR SQ MATRIX IN VECTOR MODE
! NCD     = TOTAL NUMBER OF DOF IN CONSTRAINT EQUATION
! N_D_FRE = TOTAL NUMBER OF DEGREES OF FREEDOM
! NDX(I)  = SYS DOF NOS OF CONSTRAINT TERM I
! NR      = REDUNDANT DEGREE OF FREEDOM = NDX(1)
! A(J)    = NORMALIZED COEFF OF (J+1) TERM, A0 = 1.0

!                  S*D=C SYSTEM PARTITIONED
! C = SYSTEM COLUMN MATRIX  | SII | SIR | SID | |DI| |CI|
! S = SYSTEM SQUARE MATRIX  |-----|-----|-----| |--| |--|
! I-INDEPENDENT             | SRI | SRR | SRD |*|DR|=|CR|
! R-REDUNDANT               |-----|-----|-----| |--| |--|
! D-DEPENDENT               | SDI | SDR | SDD | |DD| |CD|

  IF ( I_BUG > 0 ) PRINT *, 'Entering SKY_LCE_SYM'
  IF ( ANY (NDX < 1)       ) STOP 'SKY_LCE_SYM: BAD CONSTRAINT INDEX'
  IF ( ANY (NDX > N_D_FRE) ) STOP 'SKY_LCE_SYM: BAD CONSTRAINT INDEX'

  E     = A (NCD)            ! INITIALIZE
  NR    = NDX (1)
  SRR   = S (I_DIAG (NR) )
  SRRP1 = SRR + ONE
  CR    = C (NR)
  ESC   = E * SRRP1 - CR

!  FIND TOP AND BOTTOM OF THE SKYLINE OF THE NR COLUMN
  CALL SKY_TOP_BOT (I_DIAG, NR, ITOP, IBOT)

!  FORM MODIFIED COLUMN MATRIX, CX = CX - E*SXR
  IF ( E /= ZERO ) THEN
    DO I = ITOP, IBOT
      CALL SUB_IN_SKY (I_DIAG, I, NR, I_NR_V)
      IF ( I_NR_V /= 0 ) C (I) = C (I) - E * S (I_NR_V)
    END DO
  ENDIF

!  ADDITIONAL COLUMN CHANGES FOR CD AND CR
  IF ( NCD > 1 .AND. ESC /= ZERO ) THEN
    DO K = 2, NCD
      C (NDX (K) ) = C (NDX (K) ) + A (K - 1) * ESC
    END DO
  ENDIF
  C (NR) = E  !b add averaging data

!    *** SQUARE MATRIX COLUMN MODIFICATIONS ***
  IF ( NCD > 1)  THEN

!    FORM SID,  BEGIN SDD
    DO K = 2, NCD
      J = NDX (K)
      DO I = ITOP, IBOT
        CALL SUB_IN_SKY (I_DIAG, I, J, I_J_V)
        CALL SUB_IN_SKY (I_DIAG, I, NR, I_NR_V)
        IF ( I_NR_V /= 0 ) THEN
          IF ( S (I_NR_V) /= ZERO ) THEN
            IF ( I_J_V == 0 ) STOP 'INVALID ADDRESS, SKY_LCE_SYM'
!             DO  NOT MODIFY NR  COLUMN, MAY NEED IT LATER
            SKIP = .FALSE.
            DO LL = K + 1, NCD
              IF ( I == NDX (LL) ) SKIP = .TRUE.
            END DO
            IF ( (.NOT. SKIP) .AND. I /= NR ) &
              S (I_J_V) = S (I_J_V) - S (I_NR_V) * A (K - 1)
          END IF !  S (I_NR_V)
        END IF ! I_NR_V
      END DO
    END DO

!  COMPLETE SDD
    DO K = 2, NCD
      I = NDX (K)
      DO L = 2, K
        J = NDX (L)
        CALL SUB_IN_SKY (I_DIAG, I, J, I_J_V)
        CALL SUB_IN_SKY (I_DIAG, I, NR, I_NR_V)
        IF ( I_J_V == 0 .OR. I_NR_V == 0 )  &
             STOP 'INVALID ADDRESS, SKY_LCE_SYM'
        S (I_J_V) = S (I_J_V) + SRRP1 * A (K - 1) * A (L - 1) &
                - A (L - 1) * S (I_NR_V)
      END DO
    END DO
  ENDIF
!b XXX end LCE1 begin LCE2, split into 2 routines
!        *** INSERT CONSTRAINT EQUATION ***
! WARNING: NEXT LOOP NOT VALID FOR COUPLED CONSTRAINTS
  DO I = ITOP, IBOT
    CALL SUB_IN_SKY (I_DIAG, I, NR, I_NR_V)
    IF ( I_NR_V /= 0 ) S (I_NR_V) = ZERO
  END DO
  S (I_DIAG (NR)) = ONE   !b add averaging data
  IF ( NCD > 1 ) THEN
    DO K = 2, NCD
      I = NDX (K)
      CALL SUB_IN_SKY (I_DIAG, I, NR, I_NR_V)
      IF ( I_NR_V == 0 ) STOP 'INVALID ADDRESS, SKY_LCE_SYM'
      S (I_NR_V) = A (K - 1)  !b add averaging data
    END DO
!   *** MODIFICATIONS COMPLETED, CHECK DIAGONAL ***
    DO K = 2, NCD
      I = NDX (K)
      IF ( S (I_DIAG (I)) <= ZERO ) THEN ; WRITE (6,  * ) &
        'WARNING: NEGATIVE DIAGONAL FOR CONSTRAINT SET', NDX
        N_WARN = N_WARN + 1 ! INCREMENT WARNING
      END IF
    END DO
  ENDIF
  IF ( I_BUG > 0 ) PRINT *, 'Leaving SKY_LCE_SYM'
END SUBROUTINE SKY_LCE_SYM
                          
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
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 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 = DFREE 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 SKY_FUL_SYM (S, I_DIAG, C)
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
! COPY SKYLINE S INTO SYMMETRIC FULL S AND PRINT WITH C
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: S (N_COEFF), C (N_D_FRE)
 INTEGER,  INTENT(IN) :: I_DIAG (N_D_FRE)
!          Automatic Arrays
 REAL(DP) :: FULL (N_D_FRE, N_D_FRE)
 INTEGER  :: I, ID, IJV, ITOP, J

! N_COEFF   = NO COEFF IN SQ MATRIX = I_DIAG(N_D_FRE)
! N_D_FRE   = TOTAL NO OF DOF IN SYSTEM

  FULL = 0.d0
  DO I = 1, N_D_FRE
    DO J = I, N_D_FRE
      !b FULL (I, J) = 0.d0
!      FIND S(I, J) IN S VECTOR
      ID = MAX0 (I, J)
      IJV = I_DIAG (ID) - IABS (I - J)
      ITOP = 1
      IF ( ID > 1 ) ITOP = I_DIAG (ID-1) + 1
!      IS IT OUTSIDE SKYLINE AND THUS ZERO?
      IF ( IJV >= ITOP ) THEN
        FULL (I, J) = S (IJV)
        FULL (J, I) = FULL (I, J)
      END IF
    END DO
  END DO
  PRINT *,'SKYLINE COEFF VALUES FOLLOW'
  CALL RPRINT (FULL, N_D_FRE, N_D_FRE, 1)
  PRINT *,'COLUMN  COEFF VALUES FOLLOW'
  CALL RPRINT (C, N_D_FRE, 1, 1)
END SUBROUTINE SKY_FUL_SYM
                          
SUBROUTINE SKY_STORE_SQ_SYM (LTFREE, INDEX, I_DIAG, S, SS)
! * * * * * * * * * * * * * * * * * * * * * * * * *
!   STORE ELEMENT SQUARE MATRIX TO SYSTEM SQUARE
!    MATRIX STORED IN SYMMETRIC SKYLINE MODE
! * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 INTEGER,  INTENT(IN)    :: LTFREE
 INTEGER,  INTENT(IN)    :: INDEX (LTFREE), I_DIAG (N_D_FRE)
 REAL(DP), INTENT(IN)    :: S (LTFREE, LTFREE)
 REAL(DP), INTENT(INOUT) :: SS (N_COEFF)
 INTEGER   :: I, J, JTEMP, NDX_I, NDX_J, NDX_V
 
! N_COEFF   = NO COEFF IN SQ MATRIX = I_DIAG(N_D_FRE)
! N_D_FRE   = TOTAL NO OF DOF IN SYSTEM
! LTFREE    = NUMBER OF ELEMENT DEGREES OF FREEDOM
! INDEX(I)  = SYS DOF NO OF ELEMENT DOF I
! I_DIAG(I) = LOCATION OF DIAGONAL OF I-TH EQ
! S         = ELEMENT SQUARE MATRIX
! SS        = SYS SQ MATRIX IN SKYLINE VECTOR MODE
!
!  LOOP OVER ELEMENT COEFFICIENTS
  DO J = 1, LTFREE
    NDX_J = INDEX (J)
!    ALLOW FOR OMITTED NODES
    IF ( NDX_J > 0 ) THEN
      JTEMP = I_DIAG (NDX_J) - NDX_J
      DO I = 1, LTFREE
        NDX_I = INDEX (I)
        IF ( NDX_I <= NDX_J .AND. NDX_I > 0 ) THEN
          IF ( (I == J) .AND. (S (I, J) < 0.d0)) THEN
            PRINT *,'BAD DIAGONAL ELEMENT TERM = ', S (I, J)
            PRINT *,'ELEMENT ROW = ', I
            PRINT *,'SYSTEM  ROW = ', NDX_I
            STOP    'NEGATIVE ELEMENT DIAGONAL, SKY_STORE_SQ_SYM'
          END IF
!          FIND SYSTEM COEFF IN VECTOR S
          NDX_V = JTEMP + NDX_I
!         NDX_V = I_DIAG ( AMAX0 (NDX_I, NDX_J) )
!              - IABS (NDX_J-NDX_I)
          SS (NDX_V) = SS (NDX_V) + S (I, J)
        ENDIF
      END DO
    ENDIF
  END DO
END SUBROUTINE SKY_STORE_SQ_SYM
                               
SUBROUTINE SAVE_SKY_ROWS_GEN (I_BC, S_SKY, CC, I_DIAG, S_SKY_LOWER)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  SAVE INDEPENDENT REACTION EQUATIONS FROM A GENERAL SKYLINE ROW
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
Use System_Constants ! N_D_FRE, N_COEFF, NOT_SYM, MAX_NP, SYMMETRIC
Use Interface_Header
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: S_SKY (N_COEFF), CC (N_D_FRE)  ! always used
 REAL(DP), INTENT(IN) :: S_SKY_LOWER (NOT_SYM)          ! may be used
 INTEGER,  INTENT(IN) :: I_BC (MAX_NP), I_DIAG (N_D_FRE) ! always used
!          Automatic Arrays
 INTEGER :: INDEX (N_G_DOF)  ! equation numbers at a node
 INTEGER :: KODES (N_G_DOF)  ! bc flags at a node
 INTEGER :: J                ! node loop
 INTEGER :: IG               ! dof  loop at node
 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
                                                             
! N_COEFF = NUMBER OF TERMS IN UPPER STIFFNESS
! N_G_DOF = NUMBER OF DOF PER NODE
! N_REACT = SEQUENTIAL UNIT TO STORE REACTION DATA
! N_D_FRE = TOTAL NUMBER OF EQUATIONS
! NOT_SYM = NUMBER OF TERMS IN LOWER STIFFNESS, 0 OR N_COEFF
                                                            
  IF ( N_REACT > 0 ) THEN ; REWIND (N_REACT)
  ELSE ; STOP 'NO UNIT FOR REACTION STORAGE'
  END IF
!  LOOP OVER EQUATIONS FOR ESSENTIAL BC FLAG
  DO J = 1, MAX_NP                             ! node loop
    IF ( I_BC (J)  > 0 ) THEN                   ! bc applied
      CALL PT_CODE (J, I_BC (J), KODES) ! get kodes
      INDEX = GET_INDEX_AT_PT (J)              ! get dof at node
      DO IG = 1, N_G_DOF                       ! dof loop
        IF ( KODES (IG) < 0 ) THEN             ! component IG has bc
! XXX                   >
          ROW = INDEX (IG)                     ! extract row number
!           FIND COLUMN NUMBER BOUNDS FOR NON-ZERO COEFFICIENTS IN ROW
          COL_START = FIRST_COL_AT_SKY_ROW (ROW, I_DIAG)
          COL_STOP  =  LAST_COL_AT_SKY_ROW (ROW, I_DIAG)
!           CORRECT FOR ANY INITIALLY NULL ROWS (GRAB DIAGONAL)
          IF ( COL_START > ROW ) COL_START = ROW
          IF ( COL_STOP  < ROW ) COL_STOP  = ROW
!           WRITE NODE, PARAMETER, RANGE OF NON-ZEROS TO STORAGE
          WRITE (N_REACT) J, IG, COL_START, COL_STOP
!           SAVE ROW OF NON_ZEROS IN SYSTEM EQUATIONS
          DO COL = COL_START, COL_STOP                   ! column loop
            IN_SKY = LOCATION_IN_SKY (ROW, COL, I_DIAG)  ! locate 
              IF ( IN_SKY > 0 ) THEN                     ! non-zero ?
                IF ( ROW <= COL ) THEN                   ! in upper 
                  WRITE (N_REACT) S_SKY (IN_SKY)         ! upper value
                ELSE ! ROW > COL,                        ! in lower
                  IF ( SYMMETRIC ) THEN                  ! symmetric ?
                    WRITE (N_REACT) S_SKY (IN_SKY)       ! upper value
                  ELSE                                   ! non_sym
                    WRITE (N_REACT) S_SKY_LOWER (IN_SKY) ! lower value
                  END IF
                END IF
              ELSE
                WRITE (N_REACT) 0.d0                      ! zero term
              END IF
          END DO
          WRITE (N_REACT) CC (ROW)                       ! save row 
        END IF
      END DO
    END IF
  END DO
END SUBROUTINE SAVE_SKY_ROWS_GEN
                                
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)
Use system_constants ! N_D_FRE
 IMPLICIT NONE
 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) RESULT (COLUMN)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!    FIND LAST COLUMN IN A ROW IN A SKYLINE STORAGE
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!   (ALSO LAST ROW IN A COLUMN IN A SKYLINE STORAGE)
Use system_constants ! N_D_FRE
! Use Interface_Header
 IMPLICIT NONE
 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 FIRST_IN_SKY_ROW_OR_COL (NUMBER, I_DIAG) RESULT (IN_SKY)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! FIND VECTOR LOCATION OF FIRST NON-ZERO COEFFICIENT
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use system_constants ! N_D_FRE
Use Interface_Header
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: NUMBER            ! row or column number
 INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE)  ! skyline diagonals
 INTEGER :: IN_SKY          ! location in skyline
 INTEGER :: COL, LOCATION   ! loop
                                              
! N_D_FRE = TOTAL NUMBER OF EQUATIONS
    IN_SKY = I_DIAG(1)                   ! IF NUMBER = 1, maybe 0
    IF ( NUMBER > 1 ) THEN                        ! general case
      IN_SKY = 0                                  ! maybe null
      IF ( I_DIAG(NUMBER) > I_DIAG(NUMBER-1) ) THEN ! not null
        IN_SKY = I_DIAG(NUMBER-1) + 1              ! search over
      ELSE                                        ! zero diagonal
        DO COL = NUMBER, N_D_FRE                  ! check columns
          LOCATION = LOCATION_IN_SKY (NUMBER, COL, I_DIAG)
          IF ( LOCATION > 0 ) THEN                ! possible non-zero
            IN_SKY = LOCATION ; EXIT              ! search over
          END IF
        END DO
      END IF
    END IF
END FUNCTION FIRST_IN_SKY_ROW_OR_COL
                                    
FUNCTION LAST_IN_SKY_ROW_OR_COL (NUMBER, I_DIAG) RESULT (IN_SKY)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! FIND VECTOR LOCATION OF LAST NON-ZERO COEFFICIENT
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use system_constants ! N_D_FRE
Use Interface_Header
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: NUMBER              ! row or column number
 INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE)    ! skyline diagonals
 INTEGER :: IN_SKY            ! location in skyline
 INTEGER :: COLUMN, LOCATION  ! looping data
                                                        
! N_D_FRE = TOTAL NUMBER OF EQUATIONS
    IN_SKY = I_DIAG(NUMBER)               ! upper limit, maybe 0
    DO COLUMN = N_D_FRE, 1, -1           ! backward column loop
      LOCATION = LOCATION_IN_SKY (NUMBER, COLUMN, I_DIAG)
      IF ( LOCATION > 0 ) THEN           ! non-zero possible
        IN_SKY = LOCATION ; EXIT         ! search over
      END IF
    END DO
END FUNCTION LAST_IN_SKY_ROW_OR_COL
                                   
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
Use system_constants ! N_D_FRE
 IMPLICIT NONE
 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
                            
SUBROUTINE SKY_FUL_GEN (S, I_DIAG, C, S_LOWER)
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
! COPY SKYLINE VALUES INTO FULL FORM AND PRINT WITH C
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants ! N_D_FRE, N_COEFF, NOT_SYM, SYMMETRIC
Use Interface_Header
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: S (N_COEFF), C (N_D_FRE)
 REAL(DP), INTENT(IN) :: S_LOWER (NOT_SYM)
 INTEGER,  INTENT(IN) :: I_DIAG (N_D_FRE)
!          Automatic Arrays
 REAL(DP) :: FULL (N_D_FRE, N_D_FRE)
 INTEGER  :: I, IJV, J
                               
! N_D_FRE = TOTAL NUMBER OF EQUATIONS
  FULL = 0.
  DO I = 1, N_D_FRE                          ! row loop
    DO J = I, N_D_FRE                        ! (upper) column loop
!      FIND S(I, J) IN S VECTOR
      IJV = LOCATION_IN_SKY (I, J, I_DIAG)    ! position
!      IS IT OUTSIDE SKYLINE AND THUS ZERO?
      IF ( IJV > 0 ) THEN                    ! stored ?
        FULL (I, J) = S (IJV)                ! copy upper value
        IF ( I /= J ) THEN                   ! in lower column
          IF ( SYMMETRIC ) THEN              ! symmetric ?
            FULL (J, I) = S (IJV)            ! copy upper value
          ELSE
            FULL (J, I) = S_LOWER (IJV)      ! copy lower value
          END IF
        END IF
      END IF
    END DO
  END DO
  PRINT *,'SKYLINE COEFF VALUES FOLLOW'
  CALL RPRINT (FULL, N_D_FRE, N_D_FRE, 1)    ! print full system
  PRINT *,'COLUMN  COEFF VALUES FOLLOW'
  CALL RPRINT (C, N_D_FRE, 1, 1)
END SUBROUTINE SKY_FUL_GEN
                          
SUBROUTINE SKY_FUL_LOCATIONS (I_DIAG)
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
! PRINT SKYLINE STORAGE LOCATIONS IN FULL MATRIX FORM
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants ! N_D_FRE
Use Interface_Header
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE)
!          Automatic Arrays
 INTEGER :: FULL (N_D_FRE, N_D_FRE)
 INTEGER :: I, IJV, J
                               
! N_D_FRE = TOTAL NUMBER OF EQUATIONS
  DO I = 1, N_D_FRE                            ! row loop
    DO J = I, N_D_FRE                          ! upper col loop
!      FIND S(I, J) LOCATION IN S VECTOR(S)
      IJV = LOCATION_IN_SKY (I, J, I_DIAG)      ! may be 0
      FULL (I, J) =  IJV ; FULL (J, I) =  IJV  ! set location
      IF ( I == J ) FULL (J, J) = -IJV         ! flag diagonal
    END DO
  END DO
  PRINT *,'SKYLINE COEFF LOCATIONS FOLLOW (NEG ON DIAGONAL)'
  CALL IPRINT (FULL, N_D_FRE, N_D_FRE)         ! print locations
END SUBROUTINE SKY_FUL_LOCATIONS
                                
SUBROUTINE SHOW_SKY_FUL (I_DIAG)
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
!          SHOW SKYLINE AS A FULL CHARACTER MATRIX
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants ! N_D_FRE
Use Interface_Header
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE)
!          Automatic Arrays
 CHARACTER(LEN=1)    :: FULL (N_D_FRE, N_D_FRE)  ! work area
 CHARACTER(LEN=1)    :: TEN(10)                  ! work area
 INTEGER             :: I, I_J_V, J              ! loops
 INTEGER             :: FAKE                     ! avoid warnings
                                                                 
! N_D_FRE = TOTAL NUMBER OF EQUATIONS
  FAKE = N_D_FRE                                  ! avoid warnings
  TEN = (/ '1', '2', '3', '4', '5', '6', '7', '8', '9', ':' /)
  PRINT *, 'SKYLINE IN FULL MATRIX:'
  select case (N_D_FRE)
    case (1:9);   print *, '   ', ten(1:FAKE)
    case (10:19); print *, '   ', ten, ten(1:(FAKE-10))
    case (20:29); print *, '   ', ten, ten, ten(1:(FAKE-20))
    case (30:39); print *, '   ', ten, ten, ten, ten(1:(FAKE-30))
    case default; print *, '   ', ten, ten, ten, ten, ten
  end select
  FULL = 'o'                        ! zero symbol
  DO I = 1, N_D_FRE                 ! row loop
    DO J = I, N_D_FRE               ! upper column in row
!      FIND S(I, J) IN S VECTOR
      I_J_V = LOCATION_IN_SKY (I, J, I_DIAG)      ! may be 0
!      IS IT OUTSIDE SKYLINE AND THUS ZERO?
      IF ( I_J_V > 0 ) THEN                      ! is in skyline
        FULL (I, J) = 'U' ; FULL (J, I) = 'L'    ! upper or lower
      END IF
      IF ( I == J ) FULL (J,J) = '\\' ! single \ most systems
    END DO
  END DO
  DO I = 1, N_D_FRE                        ! row
    J = I ; IF ( J > 9 ) J = MOD(J,10)     ! row id
    PRINT *, J, ' ', FULL(I,:)                  ! id and columns in row
  END DO
END SUBROUTINE SHOW_SKY_FUL
                           
SUBROUTINE SKY_TYPE_1_GEN (LOC, VALUE, S, C, I_DIAG, S_LOWER)
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
! APPLY TYPE 1 MODIFICATION TO SYMMERTIC OR UNSYMMETRIC
! SKYLINE EQS,            S*D = C,  D(LOC) = VALUE
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants ! N_D_FRE, N_COEFF, NOT_SYM, SYMMETRIC
Use Interface_Header
 IMPLICIT NONE
 INTEGER,  INTENT(IN)    :: I_DIAG (N_D_FRE), LOC
 REAL(DP), INTENT(IN)    :: VALUE
 REAL(DP), INTENT(INOUT) :: S (N_COEFF), S_LOWER (NOT_SYM)
 REAL(DP), INTENT(INOUT) :: C (N_D_FRE)
 REAL(DP) :: AVE
 INTEGER  :: I, INV
                                           
! N_D_FRE   = DFREE NUMBER OF EQUATIONS
! N_COEFF   = NUMBER OF NON-ZERO COEFFICIENTS IN S
! NOT_SYM   = NUMBER OF NON-ZERO COEFFICIENTS IN S_LOWER
! LOC       = DOF NUMBER OF CONSTRAINED PARAMETER
! VALUE     = GIVEN VALUE OF DOF NUMBER LOC
! S         = UPPER SYSTEM MATRIX IN SKYLINE STORAGE MODE
! S_LOWER   = LOWER SYSTEM MATRIX IN SKYLINE STORAGE MODE
! C         = FULL COLUMN MATRIX
! I_DIAG     = POINTER TO DIAGONAL COEFFICIENT IN S and S_LOWER
! SYMMETRIC = true if symmetric
! AVE       = AVERAGE DIAGONAL VALUE IN CONSTRAINT
!
  IF ( LOC < 1 .OR. LOC > N_D_FRE ) THEN
    PRINT *,'WARNING: IMPOSSIBLE (FATAL) CONSTRAINT DOF = ', LOC
    PRINT *,'ASSOCIATED CONSTRAINT VALUE WAS = ', VALUE         
    N_WARN = N_WARN + 1
    STOP 'BAD LOC, SKY_TYPE_1_GEN'
  END IF
  AVE = S (I_DIAG(LOC))                         ! scaling
  if ( ave < 1.e-6) ave = 1.d0  !b
!  SUBTRACT THE LOC_TH COLUMN*VALUE FROM RHS
!  AND ZERO THE LOC_TH ROW
  DO I = 1, N_D_FRE                            ! run down the rows
!    FIND S(I, LOC) IN S VECTOR
    INV = LOCATION_IN_SKY (I, LOC, I_DIAG)
!    IS IT OUTSIDE SKYLINE AND THUS ZERO?
    IF ( INV > 0 ) THEN                  ! may not be zero, get it
      IF ( SYMMETRIC ) THEN                ! ALL IN UPPER
        C (I) = C (I) - VALUE * S (INV)    ! subtract from rhs
        S (INV) = 0.d0                      ! zero row and column
      ELSE                                 ! may be in upper or lower
        IF ( I <= LOC ) THEN                 ! use upper
          C (I) = C (I) - VALUE * S (INV)    ! subtract from rhs
          S (INV) = 0.d0                      ! zero column
          S_LOWER (INV) = 0.d0                ! zero row
        ELSE                                    ! use lower
          C (I) = C (I) - VALUE * S_LOWER (INV) ! subtract from rhs
          S_LOWER (INV) = 0.d0                   ! zero column
          S (INV) = 0.d0                         ! zero row
        END IF
      END IF
    END IF
  END DO
!  RESET THE EQUATION ROW, scaled
  IF ( .NOT. SYMMETRIC ) S_LOWER (I_DIAG(LOC)) = 0.d0
  S (I_DIAG (LOC)) = AVE
  C (LOC) = VALUE*AVE
END SUBROUTINE SKY_TYPE_1_GEN
                             
SUBROUTINE EL_HIGH (LT_INDEX, L_HIGH)
! * * * * * * * * * * * * * * * * * * * * * * * *
!  FIND SYSTEM COLUMN HEIGHTS OF AN ELEMENT
! * * * * * * * * * * * * * * * * * * * * * * * *
Use Elem_Type_Data       ! for LT_INDEX (LT_FREE)
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: LT_INDEX  (LT_FREE)
 INTEGER, INTENT(OUT) :: L_HIGH (LT_FREE)
 INTEGER :: I, MIN, NDX

! LT_FREE   = NO OF DEGREES OF FREEDOM OF ELEMENT
! LT_INDEX     = SYSTEM DOF NUMBERS OF ELEMENT PARAMETERS
! L_HIGH(I) = COLUMN HEIGHT FOR EQUATION LT_INDEX(I)

  MIN    = LT_INDEX (1)
  L_HIGH = 0
!  FIND MINIMUM LT_INDEX
  DO I = 1, LT_FREE
    NDX = LT_INDEX (I)
!    ALLOW FOR OMITTED NODES
    IF ( NDX > 0 .AND. NDX < MIN ) MIN = NDX
  END DO
!  CONVERT TO COLUMN HEIGHTS
  MIN = MIN - 1
  DO I = 1, LT_FREE
    NDX = LT_INDEX (I)
    IF ( NDX > 0 ) L_HIGH (I) = NDX - MIN
  END DO
END SUBROUTINE EL_HIGH
                     
SUBROUTINE VECTOR_MULT_SKY_GEN (S_SKY, VECTOR, I_DIAG, S_SKY_LOWER, &
                                RESULT)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! MULTIPLY A GENERAL SKYLINE ARRAY BY A VECTOR, SKY*VECTOR = RESULT
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
Use System_Constants ! N_D_FRE, N_COEFF, NOT_SYM, MAX_NP, SYMMETRIC
Use Interface_Header
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: S_SKY       (N_COEFF)     ! always used
 REAL(DP), INTENT(IN)  :: VECTOR      (N_D_FRE)     ! always used
 REAL(DP), INTENT(IN)  :: S_SKY_LOWER (NOT_SYM)     ! may be 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
! S_SKY_LOWER= LOWER HALF OF SS, IF NOT SYMMETRIC
                                                            
! 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)
    COL_STOP  =  LAST_COL_AT_SKY_ROW (ROW, I_DIAG)
!     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)  ! locate 
      IF ( IN_SKY > 0 ) THEN         ! probably non-zero 
        IF ( ROW <= COL ) THEN       ! in upper storage
          DOT = DOT + S_SKY (IN_SKY) * VECTOR (COL)        
        ELSE ! ROW > COL,            ! in lower storage
          IF ( SYMMETRIC ) THEN      ! symmetric upper value
            DOT = DOT + S_SKY (IN_SKY) * VECTOR (COL)       
          ELSE                       ! non_sym lower value
            DOT = DOT + S_SKY_LOWER (IN_SKY) * VECTOR (COL) 
          END IF
        END IF
      END IF
    END DO ! OVER COLUMNS
    RESULT (ROW) = DOT
  END DO ! OVER ROWS
END SUBROUTINE VECTOR_MULT_SKY_GEN

SUBROUTINE VECTOR_MULT_SKY_SYM (S_SKY, VECTOR, I_DIAG, RESULT)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! MULTIPLY A SYMMETRIC SKYLINE ARRAY BY A VECTOR, SKY*VECTOR = RESULT
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
Use System_Constants ! N_D_FRE, N_COEFF, NOT_SYM, MAX_NP, SYMMETRIC
Use Interface_Header
 IMPLICIT NONE
 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)
    COL_STOP  =  LAST_COL_AT_SKY_ROW (ROW, I_DIAG)
!     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)  ! 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

SUBROUTINE TYPE_1_COL_ZERO (N_RES_EQ, COEF_C_EQ, NDEX_C_EQ, CC)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!   ZERO COLUMN AT ROWS WITH EBC, FOR TRANSIENT OR DYNAMIC STEPS 
!                         OF EQS, S*D = CC
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 REAL(DP), INTENT(IN)    :: COEF_C_EQ (MAX_ACT, N_CEQ)
 INTEGER,  INTENT(IN)    :: N_RES_EQ  (MAX_ACT),             &
                            NDEX_C_EQ (MAX_ACT, N_CEQ)
 REAL(DP), INTENT(INOUT) :: CC (N_D_FRE)
 INTEGER  :: IEQ, NEQ, NTEST
                                      
! N_D_FRE         = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM
! CC              = SYS. EQ. COL. MATRIX
! N_G_DOF         = NO. PARAMETERS PER NODE
! N_RES_EQ        = NO OF CONSTRAINT EQS OF EACH TYPE
! COEF_C_EQ (I,J) = CONSTRAINT EQS COEFF I FOR EQ J
! NDEX_C_EQ (I,J) = CONSTRAINT EQS DOF NO I FOR EQ J
! N_CEQ           = NUMBER OF CONSTRAINT EQUATIONS

  IF ( I_BUG == 1 ) THEN
    PRINT *, 'Entering TYPE_1_COL_ZERO'
  END IF

!  WORK BACKWARD FROM LARGEST CONSTRAINT TYPE (DO TYPE 1 LAST)
  IEQ = SUM ( N_RES_EQ ) + 1  
    NTEST = N_RES_EQ (1) 
    IF ( NTEST > 0 ) THEN ! Apply constraint of type IC
!-->         TYPE 1     D(L1) = C1
        DO NEQ = 1, NTEST
          IEQ = IEQ - 1 
          CC (NDEX_C_EQ (1, IEQ)) = 0.d0
        END DO
    END IF
END SUBROUTINE TYPE_1_COL_ZERO
                           
