      SUBROUTINE AFACTC(MODE,NBUS)
CD
CD    THIS ROUTINE FACTORIZES AN SYMMETRIC MATRIX WITHIN THE "IBUSO",
CD    "PNODIA", "NODIA", "INDEX", "CORNO", "RELEM", "CELEM", AND
CD    "STAIRS" DATA STRUCTURE. THIS ASSUMES THAT THE GIVEN MATRIX IS
CD    STORED CONSISTENTLY WITH AN OPTIMAL ORDER OBTAINED FROM SUBROUTINE
CD    ORDER.
CD
CD    THE MEANING OF THE ARGUMENT "MODE" IS AS FOLLOWS:
CD
CD    1. MODE = 1 TOTAL FACTORIZATION REQUIRED.
CD
CD    2. MODE = 2 DETERMINE ROWS/COLUMNS TO BE 'UNREDUCED' AND
CD                AND UNREDUCE THEM.
CD
CD    3. MODE = 3 REFACTORIZE ONLY THE ROWS UNREDUCED IN (2) ABOVE
CD
CD    THE MAIN STEPS ARE GIVEN BELOW.
CD
      INCLUDE 'PARX.INS'
      INCLUDE 'CNFLAG.INS'
      INCLUDE 'CNSPC0.INS'
      INCLUDE 'CNSPDI.INS'
      INCLUDE 'CNSPRE.INS'
      INCLUDE 'CNSPC2.INS'
      INCLUDE 'CNSPVI.INS'
      INCLUDE 'CNUNIT.INS'
      INCLUDE 'CNSPWK.INS'
      COMPLEX CMPLX,PIVC,PIVR
      INTEGER*2 IORDER(NB*4),ICOLUM(NB),ICORNO(NB)
      LOGICAL LWORK(NB)
      EQUIVALENCE (RHS(1),IORDER(1)),(IORDER(NB+1),ICOLUM(1))
     &,(IORDER(NB*2+1),ICORNO(1)),(IORDER(NB*3+1),LWORK(1))
C
      LOGICAL LG
C
C
      open(IUD6,FILE='AFACTC',STATUS='OLD')
      write(IUD6,*) 'Begin Subroutine AFACTC'
      write(6,*) 'Subroutine AFACTC'
C     LDEB=(JDEB(54).NE.0)
      LDEB=.TRUE.
CD
CD    1. INITIALIZE VARIABLES DEPENDING ON VALUE OF "MODE".
CD
      DO 100 I=1,NBUS
        WORKR(I)=CMPLX(0.0,0.0)
        WORKC(I)=CMPLX(0.0,0.0)
        LWORK(I)=.FALSE.
 100  CONTINUE
C
      IBEG=1
      ILAS=NBUS-1
      INC=1
      LG=.FALSE.
C
      DO 110 I=1,NBUS
        NOROR(I)=IBUSO(I)
 110  CONTINUE
C
CD
CD    2. SET "BNO" TO BE EQUAL TO THE BUS NUMBER APPEARING LAST IN THE
CD       OPTIMAL ORDER. FOR THIS BUS NUMBER THERE IS NO CORRESPONDING
CD       ENTRY IN "RELEM", "CELEM", "CORNO", "NODIA", AND "INDEX".
CD       ONLY ITS CORRESPONDING DIAGONAL TERM APPEARS IN "STAIRS".
CD
      IBNO=IBUSO(NBUS)
C
      IF(LDEB)WRITE(IUD6,6060)(I,PNODIA(I),STAIRS(I)
     &,IBUSO(I),NODIA(I),I=1,NBUS)
CD
CD    3. START THE LOOP THAT WILL REDUCE THE "NBUS-1" ROWS/COLUMNS IN
CD       "RELEM"/"CELEM".
CD
      DO 270 II=IBEG,ILAS,INC
        IBUS=NOROR(II)
        I=PNODIA(IBUS)
        NOD=NODIA(I)
        K=INDEX(I)-1
        KA=K
        IF(LDEB)WRITE(IUD6,6040)NOD,K,II,IBUS,I
CD
CD    4. START THE LOOP WHICH UPDATES ALL THE ROWS AFFECTED BY REDUCTION
CD       OF ROW/COLUMN "I". IF ROW/COLUMN "I" CORRESPONDS TO NODE 'A',
CD       THEN THE ROWS/COLUMNS AFFECTED ARE THOSE CORRESPONDING TO THE
CD       NODES WHICH ARE CONNECTED TO NODE 'A'.
CD
        DO 250 N=1,NOD
          KA=KA+1
          NA=CORNO(KA)
          PIVC=CELEM(KA)
          PIVR=RELEM(KA)
          IF(LDEB)WRITE(IUD6,6050)KA,NA,PIVC
          IF(LDEB)WRITE(IUD6,6050)KA,NA,PIVR
CD
CD    5. UPDATE THE DIAGONAL TERM CORRESPONDING TO BUS "NA".
CD
          STAIRS(NA)=STAIRS(NA) -(PIVC)*( 1.0/STAIRS(IBUS) )*(PIVR)
CD
CD    6. CHECK IF BUS "NA" CORRESPONDS TO THE LAST BUS IN THE OPTIMAL
CD       ORDER. IF SO, THERE IS NO CORRESPONDING ROW/COLUMN IN
CD       "RELEM"/"CELEM" TO BE UPDATED.
CD
          IF(NA.EQ.IBNO)GOTO 250
          ND=PNODIA(NA)
          NODO=NODIA(ND)
          M=INDEX(ND)-1
          MA=M
CD
CD    7. START THE LOOP TO MOVE ROW/COLUMN NUMBER "ND" INTO ARRAYS
CD       "WORKR"/"WORKC" WITHIN WHICH THE UPDATING OCCURS.
CD
          DO  120 J=1,NODO
            MA=MA+1
            NR=CORNO(MA)
            WORKR(NR)=RELEM(MA)
            LWORK(NR)=.TRUE.
            WORKC(NR)=CELEM(MA)
            IF(LDEB)WRITE(IUD6,6070) NR,WORKR(NR),WORKC(NR)
  120     CONTINUE
CD
CD    8. START THE LOOP TO UPDATE ROW/COLUMN NUMBER "ND" CORRESPONDING
CD       TO BUS RESULT OF REDUCTION OF ROW NO. "I".
CD
          LA=KA
          DO 130 L= N+1,NOD
            LA=LA+1
            NR=CORNO(LA)
            LWORK(NR)=.TRUE.
            WORKR(NR)= WORKR(NR)-(PIVC)*(1.0/STAIRS(IBUS) )*(RELEM(LA))
            WORKC(NR)= WORKC(NR)-(CELEM(LA))*(1.0/STAIRS(IBUS) )*(PIVR)
            IF(LDEB)WRITE(IUD6,6080)NR,WORKR(NR),RELEM(LA),CELEM(LA)
 130      CONTINUE
CD
CD    9. PLACE UPDATED ROW/COLUMN IN "WORKR"/"WORKC" BACK INTO "CORNO",
CD       "RELEM/CELEM".
CD
C
C   TOTAL FACTORIZATION CASE
C
          JA=M
          DO 140  JB=1,NODO
            JA=JA+1
            NR=CORNO(JA)
            LWORK(NR)=.FALSE.
            NPR=PNODIA(NR)
            ICORNO(JB)=NR
            ICOLUM(JB)=NPR
            IORDER(JB)=1
 140      CONTINUE
          JJ=NODO
          LA=KA
          DO 150 JB=N+1,NOD
            LA=LA+1
            NR=CORNO(LA)
            IF(LWORK(NR))THEN
              LWORK(NR)=.FALSE.
              JJ=JJ+1
              JA=JA+1
              NPR=PNODIA(NR)
              ICORNO(JJ)=NR
              ICOLUM(JJ)=NPR
              IORDER(JJ)=1
            ENDIF
 150      CONTINUE
CD
CD    10. RESET NUMBER OF TERMS FOR ROW "ND" IN "NODIA".
CD
          NODIA(ND)=JA-M
CD
CD    11. SORT ENTRIES IN "CORNO" SO THAT BUS NUMBERS APPEAR IN THE
CD        ORDER IN WHICH THEY WERE REDUCED.
CD
          IF(LDEB)WRITE(IUD6,6130) M,JA,JJ
          IF(JJ.GT.1)THEN
            DO 160 LJ=JJ,2,-1
            IFJ=ICOLUM(LJ)
            DO 160 LI=LJ-1,1,-1
              IFI=ICOLUM(LI)
              IF(IFJ.GT.IFI)THEN
                IORDER(LJ)=IORDER(LJ)+1
              ELSE
                IORDER(LI)=IORDER(LI)+1
              ENDIF
 160        CONTINUE
            IF(LDEB)THEN
              WRITE(IUD6,6140) (ICOLUM(LL),LL=1,JJ)
              WRITE(IUD6,6150) (ICORNO(LL),LL=1,JJ)
              WRITE(IUD6,6160) (IORDER(LL),LL=1,JJ)
            ENDIF
          ENDIF
          DO 180 LI=1,JJ
            J=IORDER(LI)
            NR=ICORNO(LI)
            CORNO(M+J)=NR
            RELEM(M+J)=WORKR(NR)
            CELEM(M+J)=WORKC(NR)
            WORKR(NR)=CMPLX(0.0,0.0)
            WORKC(NR)=CMPLX(0.0,0.0)
 180      CONTINUE
C
C
         IF(LDEB)WRITE(IUD6,6090)ND,NODIA(ND)
 250  CONTINUE
 270  CONTINUE
C
      IF(LDEB)THEN
        DO 400 I=1,NBUS-1
          WRITE(IUD6,6000)I,IBUSO(I),NODIA(I)
          NE=NODIA(I)
          IS=INDEX(I)
          WRITE(IUD6,6010)(CORNO(IR),IR=IS,IS+NE-1)
          WRITE(IUD6,6020)(RELEM(IR),IR=IS,IS+NE-1)
          WRITE(IUD6,6020)(CELEM(IR),IR=IS,IS+NE-1)
 400    CONTINUE
        WRITE(IUD6,6110)
        WRITE(IUD6,6120) (STAIRS(I),I=1,NBUS)
      ENDIF
      write(IUD6,*) 'end Subroutine AFACTC'
      CLOSE(IUD6)
      RETURN
C
 6000 FORMAT(//,1X,'ROW =',I5,'OLD BUS=',I5,'NO. ELEMENTS=',I5,/)
 6010 FORMAT(1X,10I10)
 6020 FORMAT(1X,10G12.4)
 6040 FORMAT(1X,'NOD=',I5,'  K=',I5,' II=',I5,' IBUS=',I5,' I=',I5)
 6050 FORMAT(/,1X,'KA=',I5,'  NA=',I5,'  PIVC=',4G12.5)
 6051 FORMAT(/,1X,'KA=',I5,'  NA=',I5,'  PIVR=',4G12.5)
 6060 FORMAT(/,' I=',I5,' PNODIA(I)=',I5,' STAIRS(I)=',2G12.5
     &,' IBUSO(I)=',I5,' NODIA(I)=',I5)
 6070 FORMAT(1X,'1020 NR=',I5,'  WORKR=',8G12.5)
 6080 FORMAT(1X,'1030 NR=',I5,'  WORKR=',6G9.3)
 6090 FORMAT(/,1X,'ND=',I5,'  NODIA(ND)=',I5)
 6110 FORMAT(/,1X,'STAIRS')
 6120 FORMAT(/,1X,10G12.5)
 6130 FORMAT(1X,' M=',I5,' JA=',I5,' JJ=',I5)
 6140 FORMAT(1X,'ICOLUM ARRAY',/,(20I5))
 6150 FORMAT(1X,'ICORNO ARRAY',/,(20I5))
 6160 FORMAT(1X,'IORDER ARRAY',/,(20I5))
      END
