IMPLICIT NONE
INTEGER :: I, J, K, COUNT = 0, COUNT_HIER = 0
INTEGER :: MIN_DF_COLUMN, SUM_DF_COLUMN 
INTEGER :: N, M, Y, Z, ELEMENTS
CHARACTER (LEN = 1), DIMENSION (:,:), ALLOCATABLE :: FUN, RAR_FUN_ROW, &
   RAR_FUN_COL
INTEGER, DIMENSION (:), ALLOCATABLE :: DF_COLUMN, VARIABLES
INTEGER, DIMENSION (:), ALLOCATABLE :: DF_ROW, HIER, EQNS
INTEGER, DIMENSION (:), ALLOCATABLE :: ROW_FREQ
INTEGER :: ROW_FREQ_MIN, HIGHEST_HIERARCHY, ROWS_ELIMINATED 
INTEGER :: ROWNUM_MIN_FREQ, TOP_OF_HIER
CHARACTER (LEN = 30) :: MATRIX_FORMAT, COUNT_FORMAT, COL_DF_T_FORMAT, HEADER
CHARACTER (LEN = 30) :: FINAL_FUN_FORMAT
INTEGER :: COUNT_ROW_ORDER, DIAG_ELEMENT, COUNT_SUB_GROUP
INTEGER, DIMENSION (:), ALLOCATABLE :: FINAL_ROW_ORDER, FINAL_EQNS, SUB_GROUP
CHARACTER (LEN = 1), DIMENSION (:,:), ALLOCATABLE :: FINAL_FUN_MATRIX 
CHARACTER (LEN = 50) :: HIERARCHY_FORMAT, SUB_GROUP_FORMAT, READ_FORMAT
CHARACTER (LEN = 1) :: LABEL
INTEGER :: HIERARCHY_ROW_ALREADY_ELIMINATED = 0

READ_FORMAT = "( I2, 1X, I2, 1X, A)" 

OPEN (UNIT = 1, FILE = 'functionality.matrix', STATUS = 'OLD')

READ(1,*) N        
READ(1,*) M
READ(1,*) ELEMENTS

ALLOCATE (FUN (N,M), RAR_FUN_ROW(N,M), RAR_FUN_COL(N,M))
ALLOCATE (DF_COLUMN(M), VARIABLES(M))
ALLOCATE (DF_ROW(N), HIER(N), EQNS(N), ROW_FREQ(N))
ALLOCATE (FINAL_ROW_ORDER(N), FINAL_EQNS(N), SUB_GROUP(N))
ALLOCATE (FINAL_FUN_MATRIX(N,M))

DO I = 1,N
  DO J = 1,M
    FUN(I,J) = '-'
    RAR_FUN_ROW(I,J) = '-'
    RAR_FUN_COL(I,J) = '-'
  END DO
END DO

DO I = 1,ELEMENTS
  READ(1,READ_FORMAT) Y, Z, LABEL
  FUN(Y,Z) = LABEL
END DO

CLOSE(1)


DO I = 1,N
  EQNS(I) = I
END DO
DO I = 1,M
  VARIABLES(I) = I
END DO

HEADER = "(1X, /, 4X, 100(1X,I2)/)" 
FINAL_FUN_FORMAT = "(1X, I2, 1X, 100(2X,A))"
MATRIX_FORMAT = "(1X, 4(2X,A))"
SUB_GROUP_FORMAT = "(1X, 'SUB GROUP', I3, ' = ROW', I3)" 
HIERARCHY_FORMAT = "(1X, 'HIERARCHY', I3, ' = ROW', I3)"

WRITE(*,HEADER) VARIABLES
DO I = 1,N 
 WRITE(*,FINAL_FUN_FORMAT) EQNS(I), (FUN(I,J), J = 1,M)
END DO
! search for column degree of freedom of 1

CALL CALC_COLUMN_DF(N, M, DF_COLUMN, FUN, MIN_DF_COLUMN, SUM_DF_COLUMN)

DO WHILE(SUM_DF_COLUMN /= 0)
  DO WHILE(MIN_DF_COLUMN == 1)

    CALL ONE_DF_COLUMN_ARRANGE(N, M, COUNT, DF_COLUMN, FUN, RAR_FUN_ROW, &
      EQNS)

    CALL CALC_COLUMN_DF(N, M, DF_COLUMN, FUN, MIN_DF_COLUMN, SUM_DF_COLUMN)

  END DO

! search for degree of freedom greater than 1

  CALL MANY_DF_COLUMN_ARRANGE(N, M, COUNT, DF_COLUMN, MIN_DF_COLUMN, &
   FUN, RAR_FUN_ROW, COUNT_HIER, HIER, EQNS)

  CALL CALC_COLUMN_DF(N, M, DF_COLUMN, FUN, MIN_DF_COLUMN, SUM_DF_COLUMN)

! return to search for 1 degree of freedom

END DO

COL_DF_T_FORMAT = "(1X, A, I3)"
COUNT_FORMAT = "(1X, A, I4)"

! ok, now that the rows have been eliminated, lets eliminate the columns

WRITE(*,HEADER) VARIABLES
DO I = 1,N
 WRITE(*,FINAL_FUN_FORMAT) EQNS(I), (RAR_FUN_ROW(I,J), J = 1,M)
END DO

COUNT = 0
COUNT_ROW_ORDER  = 0
COUNT_SUB_GROUP = 0
DIAG_ELEMENT = 1
HIGHEST_HIERARCHY = COUNT_HIER
TOP_OF_HIER = 1
IF (COUNT_HIER == 0) HIGHEST_HIERARCHY = 1

! if the hier(element) is zero, then set it to the number of rows, because
! so that if no hierarchy exists, the algorithm will treat the entire 
! system of equations as the current hierarchy
! if there are no hierarchy pointers, then highest_hierarchy set to 1
! because hier(highest_hierarchy) would be hier(0) which doesn't exist 

DO I = 1,N
  IF (HIER(I) == 0) HIER(I) = N
END DO

IF (COUNT_HIER /= 0) THEN
  DO I = COUNT_HIER, 1, -1
    WRITE(*,HIERARCHY_FORMAT) I, EQNS(HIER(I))
  END DO
ELSE
  WRITE(*,*) 'NO HIERARCHY POINTERS' 
END IF

! establish where the higest hierarchy is or if there are any hierarchies

CALL CHECK_HIER_STATUS(N, M, RAR_FUN_ROW, HIER, COUNT_HIER, &
     TOP_OF_HIER, HIGHEST_HIERARCHY, HIERARCHY_ROW_ALREADY_ELIMINATED)

! set up a counter loop to end execution after all rows eliminated

ROWS_ELIMINATED = 0

DO WHILE (ROWS_ELIMINATED < N)

! calculate the row frequencies on the highest hierarchy

  CALL CALC_ROW_FREQ(N, M, RAR_FUN_ROW, ROW_FREQ, ROW_FREQ_MIN, & 
    ROWNUM_MIN_FREQ, TOP_OF_HIER, HIER, HIGHEST_HIERARCHY)

  CALL CHECK_HIER_STATUS(N, M, RAR_FUN_ROW, HIER, COUNT_HIER, &
     TOP_OF_HIER, HIGHEST_HIERARCHY, HIERARCHY_ROW_ALREADY_ELIMINATED)

! delete one row if the minimum frequency on the current hierarchy is 1

  IF (ROW_FREQ_MIN == 1) THEN

    CALL ELIMINATE_ONE_ROW(N, M, RAR_FUN_ROW, RAR_FUN_COL, &
      ROWNUM_MIN_FREQ, COUNT, VARIABLES, COUNT_ROW_ORDER, &
      FINAL_ROW_ORDER, DIAG_ELEMENT, HIER, HIGHEST_HIERARCHY, &
      HIERARCHY_ROW_ALREADY_ELIMINATED, TOP_OF_HIER)

  END IF

! however, if the minumum frequency on the hierarchy is greater than one
! delete the top row

  IF (ROW_FREQ_MIN /= 1 .AND. ROW_FREQ_MIN /= 0) THEN

    CALL ELIMINATE_TOP_ROW(N, M, RAR_FUN_ROW, RAR_FUN_COL, ROW_FREQ, COUNT, &
     VARIABLES, COUNT_ROW_ORDER, FINAL_ROW_ORDER, DIAG_ELEMENT, &
     COUNT_SUB_GROUP, SUB_GROUP, EQNS, TOP_OF_HIER, HIER, &
     HIGHEST_HIERARCHY, HIERARCHY_ROW_ALREADY_ELIMINATED)  

  END IF

  ROWS_ELIMINATED = ROWS_ELIMINATED + 1

END DO

! form the final rearranged matrix accounting for shuffled rows

CALL RAR_FINAL_ROW_ORDER(N, M, FINAL_ROW_ORDER, RAR_FUN_COL, &
  FINAL_FUN_MATRIX, FINAL_EQNS, EQNS)
WRITE(*,HEADER) VARIABLES
DO I = 1,N
  WRITE(*,FINAL_FUN_FORMAT) FINAL_EQNS (I), (FINAL_FUN_MATRIX(I, J), J=1,M)
END DO

! print out the relevant subgroups discoverd

DO I = 1,COUNT_SUB_GROUP
  WRITE(*,SUB_GROUP_FORMAT) I, SUB_GROUP(I)
END DO  
    
END


SUBROUTINE CALC_COLUMN_DF(N, M, DF_COLUMN, FUN, MIN_DF_COLUMN, SUM_DF_COLUMN)
IMPLICIT NONE
INTEGER :: I, J
INTEGER :: N, M
INTEGER, DIMENSION (M) :: DF_COLUMN
CHARACTER (LEN = 1), DIMENSION (N,M) :: FUN
INTEGER :: MIN_DF_COLUMN, SUM_DF_COLUMN

! calculate the degrees of freedom for each row

  DF_COLUMN(1:M) = 0

  DO I = 1,M
    DO J = 1,N
      IF (FUN(J,I) /= '-') DF_COLUMN(I) = DF_COLUMN(I) + 1
    END DO
  END DO
  
! calculate the minimum degrees of freedom and the sum of all the column
! degrees of freedom so the routine knows when to stop eliminating rows, 
! i.e. the sum is zero
! note that 1000 is an arbitrarily large number

  MIN_DF_COLUMN = 1000

  DO I = 1,M
    IF (DF_COLUMN(I) < MIN_DF_COLUMN .AND. DF_COLUMN(I) /= 0)  &
       MIN_DF_COLUMN = DF_COLUMN(I)
  END DO
        
  SUM_DF_COLUMN = SUM(DF_COLUMN)


END




SUBROUTINE ONE_DF_COLUMN_ARRANGE(N, M, COUNT, DF_COLUMN, FUN, & 
   RAR_FUN_ROW, EQNS) 
IMPLICIT NONE
INTEGER :: I, J, COUNT
INTEGER :: N, M
INTEGER, DIMENSION (M) :: DF_COLUMN
CHARACTER (LEN = 1), DIMENSION (N,M) :: FUN, RAR_FUN_ROW
INTEGER, DIMENSION (N) :: EQNS

! eliminate one row because the minimum column degree of freedom is 1

  DO I = 1,M
    IF(DF_COLUMN(I) == 1) THEN
      COUNT = COUNT + 1
      DO J = 1,N
        IF(FUN(J,I) /= '-') THEN
           RAR_FUN_ROW(N+1-COUNT,1:M) = FUN(J,1:M)
           FUN(J,1:M) = '-'
           EQNS(N+1-COUNT) = J
        END IF
      END DO
      EXIT    
    END IF  
  END DO

END


SUBROUTINE MANY_DF_COLUMN_ARRANGE(N, M, COUNT, DF_COLUMN, MIN_DF_COLUMN, & 
   FUN, RAR_FUN_ROW, COUNT_HIER, HIER, EQNS)
IMPLICIT NONE
INTEGER :: I, J, COUNT, MIN_DF_COLUMN, COUNT_HIER
INTEGER :: HIERARCHY
INTEGER :: N, M
CHARACTER (LEN = 1), DIMENSION (N,M) :: FUN, RAR_FUN_ROW
INTEGER, DIMENSION (M) :: DF_COLUMN
INTEGER, DIMENSION (N) :: HIER, EQNS

! eliminate all rows of one column and assign a hierarchy pointer because
! the minimum degree of freedom is greater than one

  HIERARCHY = 0
  DO I = 1,M   
    IF(DF_COLUMN(I) == MIN_DF_COLUMN) THEN
       DO J = 1,N
         IF(FUN(J,I) /= '-') THEN
           HIERARCHY = HIERARCHY + 1
           COUNT = COUNT + 1
           RAR_FUN_ROW(N+1-COUNT,1:M) = FUN(J,1:M)
           FUN(J,1:M) = '-'
           EQNS(N+1-COUNT) = J
             IF (HIERARCHY == 1) THEN
                COUNT_HIER = COUNT_HIER + 1
                HIER(COUNT_HIER) = N+1-COUNT
             END IF
         END IF
       END DO   
       EXIT
    END IF
  END DO

END

SUBROUTINE CALC_ROW_FREQ(N, M, RAR_FUN_ROW, ROW_FREQ, ROW_FREQ_MIN, &
     ROWNUM_MIN_FREQ, TOP_OF_HIER, HIER, HIGHEST_HIERARCHY)
IMPLICIT NONE
INTEGER :: I, J, ROW_FREQ_MIN, ROWNUM_MIN_FREQ
INTEGER :: N, M, TOP_OF_HIER, HIGHEST_HIERARCHY, ken
INTEGER, DIMENSION (N) :: ROW_FREQ, HIER
CHARACTER (LEN = 1), DIMENSION (N,M) :: RAR_FUN_ROW

! calculate the row degree of freedom in the uppermost hierarchy 

ROW_FREQ(TOP_OF_HIER:HIER(HIGHEST_HIERARCHY)) = 0

DO I = TOP_OF_HIER,HIER(HIGHEST_HIERARCHY)
  DO J = 1,M
     IF (RAR_FUN_ROW(I,J) /= '-') ROW_FREQ(I) = ROW_FREQ(I) + 1
  END DO
END DO

! find out what row has the minimum degrees of freedom
! note that 1000 is simply an arbitrarily large number

ROW_FREQ_MIN = 1000

DO I = TOP_OF_HIER,HIER(HIGHEST_HIERARCHY)
  IF (ROW_FREQ(I) /= 0 .AND. ROW_FREQ(I) < ROW_FREQ_MIN) THEN
    ROW_FREQ_MIN = ROW_FREQ(I)
    ROWNUM_MIN_FREQ = I
  END IF
END DO

END SUBROUTINE CALC_ROW_FREQ



SUBROUTINE ELIMINATE_TOP_ROW(N, M, RAR_FUN_ROW, RAR_FUN_COL, ROW_FREQ, &
   COUNT, VARIABLES, COUNT_ROW_ORDER, FINAL_ROW_ORDER, DIAG_ELEMENT, &
   COUNT_SUB_GROUP, SUB_GROUP, EQNS, TOP_OF_HIER, HIER, HIGHEST_HIERARCHY,&
   HIERARCHY_ROW_ALREADY_ELIMINATED)
IMPLICIT NONE
INTEGER :: N, M
INTEGER :: I, J, COUNT_ELIM_ROW, COUNT, L, TOP_OF_HIER, HIGHEST_HIERARCHY
INTEGER, DIMENSION (N) :: ROW_FREQ
CHARACTER (LEN = 1), DIMENSION (N,M) :: RAR_FUN_ROW, RAR_FUN_COL
INTEGER, DIMENSION (M) :: VARIABLES
INTEGER :: COUNT_ROW_ORDER, DIAG_ELEMENT, COUNT_SUB_GROUP, ROW_SHUFFLED
INTEGER, DIMENSION (N) :: FINAL_ROW_ORDER, SUB_GROUP, EQNS, HIER
INTEGER :: FREQ_HIERARCHY_ROW, HIERARCHY_ROW_ALREADY_ELIMINATED

! eliminate all the columns with entries in the top row because the row 
! because the minimum row degree of freedom is greater than 1

COUNT_ELIM_ROW = 0

DO I = TOP_OF_HIER, HIER(HIGHEST_HIERARCHY)
  ROW_SHUFFLED = 0
  DO L = TOP_OF_HIER, HIER(HIGHEST_HIERARCHY)
    IF (FINAL_ROW_ORDER(L) == I) ROW_SHUFFLED = 1
  END DO
  IF (ROW_SHUFFLED /= 1) THEN
    IF (ROW_FREQ(I) /= 0) COUNT_ELIM_ROW = COUNT_ELIM_ROW + 1
    IF (COUNT_ELIM_ROW == 1) THEN
      DO J = 1,M
         IF (RAR_FUN_ROW(I,J) /= '-') THEN
           COUNT = COUNT + 1
           RAR_FUN_COL(1:N, COUNT) = RAR_FUN_ROW(1:N, J)
           RAR_FUN_ROW(1:N, J) = '-'
           VARIABLES(COUNT) = J
         END IF
      END DO
      COUNT_ROW_ORDER = COUNT_ROW_ORDER + 1
      FINAL_ROW_ORDER(COUNT_ROW_ORDER) = I
      DIAG_ELEMENT = I
      COUNT_SUB_GROUP = COUNT_SUB_GROUP + 1
      SUB_GROUP(COUNT_SUB_GROUP) = EQNS(I)
    END IF
  END IF

END DO

! account for the possibility that eliminating one row above the 
! hierarchy pointer may simultaneously eliminate the hierarchy row
! but not if it has already been eliminated because it was the minimum
! frequency row with a frequency of one

IF (HIERARCHY_ROW_ALREADY_ELIMINATED == 0) THEN

FREQ_HIERARCHY_ROW = 0

DO I = 1,M
  IF (RAR_FUN_ROW(HIER(HIGHEST_HIERARCHY), I) /= '-') THEN
    FREQ_HIERARCHY_ROW = FREQ_HIERARCHY_ROW + 1
  END IF
END DO
        
IF (FREQ_HIERARCHY_ROW == 0) THEN
  COUNT_ROW_ORDER = COUNT_ROW_ORDER + 1
  FINAL_ROW_ORDER(COUNT_ROW_ORDER) = HIER(HIGHEST_HIERARCHY)
END IF

END IF

END SUBROUTINE ELIMINATE_TOP_ROW


     
SUBROUTINE ELIMINATE_ONE_ROW(N, M, RAR_FUN_ROW, RAR_FUN_COL, &
  ROWNUM_MIN_FREQ, COUNT, VARIABLES, COUNT_ROW_ORDER, FINAL_ROW_ORDER, &
  DIAG_ELEMENT, HIER, HIGHEST_HIERARCHY, HIERARCHY_ROW_ALREADY_ELIMINATED, &
  TOP_OF_HIER) 
IMPLICIT NONE
INTEGER :: N, M
INTEGER :: I, J, COUNT, COUNT_ELIM_ROW, ROWNUM_MIN_FREQ
CHARACTER (LEN = 1), DIMENSION (N,M) :: RAR_FUN_ROW, RAR_FUN_COL
INTEGER, DIMENSION (M) :: VARIABLES
INTEGER :: COUNT_ROW_ORDER, DIAG_ELEMENT, HIGHEST_HIERARCHY
INTEGER :: FREQ_HIERARCHY_ROW
INTEGER, DIMENSION (N) :: FINAL_ROW_ORDER, HIER
INTEGER :: HIERARCHY_ROW_ALREADY_ELIMINATED
INTEGER, DIMENSION (N) :: PRE_ROW_DF, POST_ROW_DF, POST_MINUS_PRE
INTEGER :: TOP_OF_HIER

! check for the condition that eliminating one row may eliminate 
! another.  The possiblilty of eliminating the hierarchy pointer
! is checked for later

PRE_ROW_DF(TOP_OF_HIER:HIER(HIGHEST_HIERARCHY)) = 0
DO I = TOP_OF_HIER, (HIER(HIGHEST_HIERARCHY)-1)
  DO J = 1,M
    IF (RAR_FUN_ROW(I,J) /= '-') PRE_ROW_DF(I) = PRE_ROW_DF(I) + 1
  END DO
END DO

COUNT_ELIM_ROW = 0

! eliminate one row because the minimum row degree of freedom is 1 

DO I = 1, M
  IF (RAR_FUN_ROW(ROWNUM_MIN_FREQ, I) /= '-') THEN
     COUNT_ELIM_ROW = COUNT_ELIM_ROW + 1
     IF (COUNT_ELIM_ROW == 1) THEN
        COUNT = COUNT + 1
        RAR_FUN_COL(1:N, COUNT) = RAR_FUN_ROW(1:N, I)
        RAR_FUN_ROW(1:N, I) = '-'
        VARIABLES(COUNT) = I
        COUNT_ROW_ORDER = COUNT_ROW_ORDER + 1

! shuffle a row if it is not on the current diagonal, yet above the current
! hierarchy pointer

        IF (ROWNUM_MIN_FREQ >= DIAG_ELEMENT) THEN
          FINAL_ROW_ORDER(COUNT_ROW_ORDER) = ROWNUM_MIN_FREQ
          DIAG_ELEMENT = DIAG_ELEMENT + 1
        ELSE
          FINAL_ROW_ORDER(COUNT_ROW_ORDER) = I
        ENDIF

     END IF
  END IF
END DO

POST_ROW_DF(TOP_OF_HIER:HIER(HIGHEST_HIERARCHY)) = 0
DO I = TOP_OF_HIER, (HIER(HIGHEST_HIERARCHY)-1)
  DO J = 1,M
    IF (RAR_FUN_ROW(I,J) /= '-') POST_ROW_DF(I) = POST_ROW_DF(I) + 1
  END DO
END DO

POST_MINUS_PRE = POST_ROW_DF - PRE_ROW_DF
DO I = TOP_OF_HIER, (HIER(HIGHEST_HIERARCHY)-1)
IF (POST_MINUS_PRE(I) < 0 .AND. I /= ROWNUM_MIN_FREQ) THEN
  IF (POST_ROW_DF(I) == 0) THEN
    COUNT_ROW_ORDER = COUNT_ROW_ORDER + 1
    FINAL_ROW_ORDER(COUNT_ROW_ORDER) = I
  END IF
END IF
END DO

! account for the possibility that eliminating one row above the
! hierarchy pointer may simultaneously eliminate the hierarchy row
! but not if it has already been eliminated because it was the minimum
! frequency row with a frequency of one

IF (ROWNUM_MIN_FREQ == HIER(HIGHEST_HIERARCHY)) THEN
  HIERARCHY_ROW_ALREADY_ELIMINATED = 1
END IF

IF (HIERARCHY_ROW_ALREADY_ELIMINATED == 0) THEN

FREQ_HIERARCHY_ROW = 0

DO I = 1,M
  IF (RAR_FUN_ROW(HIER(HIGHEST_HIERARCHY), I) /= '-') THEN
    FREQ_HIERARCHY_ROW = FREQ_HIERARCHY_ROW + 1
  END IF
END DO

IF (FREQ_HIERARCHY_ROW == 0) THEN
  COUNT_ROW_ORDER = COUNT_ROW_ORDER + 1
  FINAL_ROW_ORDER(COUNT_ROW_ORDER) = HIER(HIGHEST_HIERARCHY)
END IF

END IF

END SUBROUTINE ELIMINATE_ONE_ROW


SUBROUTINE CHECK_HIER_STATUS(N, M, RAR_FUN_ROW, HIER, COUNT_HIER, &
  TOP_OF_HIER, HIGHEST_HIERARCHY, HIERARCHY_ROW_ALREADY_ELIMINATED) 
IMPLICIT NONE 
INTEGER :: N, M
INTEGER :: I, J, HIGHEST_HIERARCHY
CHARACTER (LEN = 1), DIMENSION (N,M) :: RAR_FUN_ROW
INTEGER, DIMENSION (N) :: HIER  
INTEGER :: COUNT_HIER, TOP_OF_HIER, HIER_STATUS
INTEGER :: HIERARCHY_ROW_ALREADY_ELIMINATED

! determine which hierarchy pointer is currently the highest and
! if there are no hierarchy pointers set the current hierarchy
! to be the last row (entire matrix)

HIER_STATUS = 0
DO I = TOP_OF_HIER, HIER(HIGHEST_HIERARCHY)
  DO J = 1, M
    IF (RAR_FUN_ROW(I, J) /= '-') HIER_STATUS = 1
  END DO
END DO

IF (HIER_STATUS == 0) THEN
  TOP_OF_HIER = HIER(HIGHEST_HIERARCHY) + 1
  COUNT_HIER = COUNT_HIER - 1
  HIGHEST_HIERARCHY = COUNT_HIER
  IF (COUNT_HIER == 0) HIGHEST_HIERARCHY = N
  HIERARCHY_ROW_ALREADY_ELIMINATED = 0
END IF

END SUBROUTINE CHECK_HIER_STATUS


SUBROUTINE RAR_FINAL_ROW_ORDER(N, M, FINAL_ROW_ORDER, RAR_FUN_COL, &
  FINAL_FUN_MATRIX, FINAL_EQNS, EQNS)
IMPLICIT NONE
INTEGER :: N, M, I, NEXT_ROW
CHARACTER (LEN = 1), DIMENSION (N,M) :: RAR_FUN_COL, FINAL_FUN_MATRIX
INTEGER, DIMENSION (N) :: FINAL_ROW_ORDER, EQNS, FINAL_EQNS

! properly order the rows in the output matrix since they may have been
! shuffled during the variable ordering algorithm

DO I = 1,N
  NEXT_ROW = FINAL_ROW_ORDER(I)
  FINAL_FUN_MATRIX(I, 1:M) = RAR_FUN_COL(NEXT_ROW, 1:M)
  FINAL_EQNS(I) = EQNS(NEXT_ROW)
END DO

END SUBROUTINE RAR_FINAL_ROW_ORDER

