      PROGRAM TWOLINK
      IMPLICIT NONE

C     Calculates the links of radius 2 of the triangle presentations (of
C     index 3 subgroups of the cases with torsion) given as
C     input. Outputs the adjacency matrices of the two-links in a form
C     that can be inputted to Magma.

C     There is 45 labels in the triangles. We denote them as follows:
C     incoming 1 = 1
C     outgoing 1 = 2
C     incoming 2 = 3
C     outgoing 2 = 4
C     incoming 3 = 5 etc.

C     The 2-link is calculated around that vertex which has incoming
C     edge with the label 1. In that vertex, the incoming edges have
C     labels 3k-2 and outgoing edges of with labels 3k. The 2-link is
C     calculated starting with the neighbour along incoming 1, then
C     outgoing 3, incoming 4, outgoin 6 and so on. We calculate the
C     1-links at the neighbouring vertices, then remove and add certain
C     edges in/between the 1-links to obtain 2-link. The resulta is a
C     graph with 780 vertices, all of index 3 (30 x part of 1-link +
C     additional edges between these).

      INTEGER I,J,L,R,M, TRIPS(45*7159,3)
      INTEGER TRIP(45,3), K, NB, COL, ROW,OK, MAX
      INTEGER KK, JJ, II, YI, XJ, TRIPNB, EDGE, CMAX
      INTEGER REMLIST(4*30), PREV(3,2), ADDEDGE(3,2,2)
      INTEGER FROW, FCOL,REMNB, LABEL, EDGESLOT
C     Adjacency matrix for the link of radius 1:
      INTEGER ADJONEA(90,90), ADJONEB(90,90), ADJONEC(90,90)
      INTEGER FDATA(30*30*4,2), LDATA(2700)
C The list of cases that we want to consider
      INTEGER UNIQS(168)
C POISTA LOPUKSI
      INTEGER ADJONE(30,30)

C     Adjacency matrix for the link of radius 2:
      INTEGER ADJTWO(30*90,30*90)
C     The same without the extra vertices, ie. final OUTPUT
      INTEGER FINAL(30*26, 30*26)

      MAX=7159
      CMAX=168

C     Open the data file on unit 3
 20   OPEN (UNIT=3, FILE='inx3sgs.txt')
      OPEN (UNIT=7, FILE='2links.txt')
      OPEN (UNIT=4, FILE='classes.txt')

C     Read the data into TRIPS & UNIQS
 25   READ (3,*) ((TRIPS(I,J), J=1,3), I=1,45*MAX)
      CLOSE(UNIT=3)
      READ(4,*) (UNIQS(I),I=1,CMAX)
      CLOSE(UNIT=4)

      WRITE(*,*) 'data reading ok'

      TRIPNB=-1

 100  TRIPNB=TRIPNB+1
      REMNB=0

      IF(TRIPNB .LT. MAX) THEN
         DO 120 I=1,CMAX
            IF(TRIPNB+1 .EQ. UNIQS(I)) THEN
               GOTO 130
            END IF
 120     CONTINUE
         GOTO 100
      ELSE
         GOTO 9999
      END IF
 

C     Read one of the groups into TRIP
 130   DO 140 I=1,45
         TRIP(I,1)=TRIPS(45*TRIPNB+I,1)
         TRIP(I,2)=TRIPS(45*TRIPNB+I,2)
         TRIP(I,3)=TRIPS(45*TRIPNB+I,3)
 140   CONTINUE

      WRITE(*,*) TRIPNB+1

C     Write the adjacency matrix for the one-link into ADJONEA, ADJONEB
C     and ADJONEC

C To ADJONEA arrives edges 3i-2, leaves 3i
C To ADJONEB arrives edges 3i, leaves 3i-1
C To ADJONEC arrives edges 3i-1, leaves 3i-2

 145   DO 148 I=1,90
 146      DO 147 J=1,90
            ADJONEA(I,J)=0
            ADJONEB(I,J)=0
            ADJONEC(I,J)=0
 147      CONTINUE
 148   CONTINUE


 150   DO 160 I=1,45
         DO 155 J=1,2
            II=TRIP(I,J)
            JJ=TRIP(I,J+1)
            YI=2*II-1
            XJ=2*JJ
            IF (MOD(II,3).EQ.1) THEN
               ADJONEA(YI,XJ)=1
               ADJONEA(XJ,YI)=1
            ELSEIF(MOD(II,3).EQ.0) THEN
               ADJONEB(YI,XJ)=1
               ADJONEB(XJ,YI)=1
            ELSE
               ADJONEC(YI,XJ)=1
               ADJONEC(XJ,YI)=1
            END IF
 155      CONTINUE
            II=TRIP(I,3)
            JJ=TRIP(I,1)
            YI=2*II-1
            XJ=2*JJ
            IF (MOD(II,3).EQ.1) THEN
               ADJONEA(YI,XJ)=1
               ADJONEA(XJ,YI)=1
            ELSEIF(MOD(II,3).EQ.0) THEN
               ADJONEB(YI,XJ)=1
               ADJONEB(XJ,YI)=1
            ELSE
               ADJONEC(YI,XJ)=1
               ADJONEC(XJ,YI)=1
            END IF
 160   CONTINUE


C       WRITE(*,*) '1-link A:'
C       WRITE(*,'(90I1)') ((ADJONEA(I,J),J=1,90),I=1,90)
C       WRITE(*,*) '1-link B:'
C       WRITE(*,'(90I1)') ((ADJONEB(I,J),J=1,90),I=1,90)
C       WRITE(*,*) '1-link C:'
C       WRITE(*,'(90I1)') ((ADJONEC(I,J),J=1,90),I=1,90)

C     We construct 2-link around a vertex of the type A


C     Take 15 copies of the 1-link adjacency matrix ADJONEB and 15
C     copies of ADJONEC on the diagonal of the two-link adjagency matrix

       DO 180 I=1,30*90
          DO 170 J=1,30*90
             ADJTWO(I,J)=0
 170      CONTINUE
 180   CONTINUE


 200  DO 300 I=1,30
         DO 250 J=1,90
            DO 240 K=1,90
               IF(MOD(I,2).EQ. 1) THEN
                  ADJTWO(90*(I-1)+K,90*(I-1)+J)=ADJONEC(K,J)
               ELSE 
                  ADJTWO(90*(I-1)+K,90*(I-1)+J)=ADJONEB(K,J)
               END IF
 240        CONTINUE
 250     CONTINUE 
 300  CONTINUE


C Remove certain edges from the graph
      
 1000 DO 1500 LABEL=1,15
 1001    DO 1499 II=1,2
C     For even i, this considers the x_i=NB  in the 1-link
C     For odd i, this considers the y_i=NB in the 1-link
            NB=2*(LABEL-1)+II
            IF(II.EQ.1) THEN
C     3*LABEL-2 into A
               EDGE=3*LABEL-2
C     out from C
               EDGESLOT=2*EDGE
            ELSE
C     3*LABEL out from A
               EDGE=3*LABEL
C     into B
               EDGESLOT=2*EDGE-1
            END IF
C     FIND all j such that there is triple (j,EDGE,k) in case of x_i
C     and triple (EDGE,j,k) in case of y_i.
C     List the triplets to PREV.

C            WRITE(*,*) EDGE, EDGESLOT

            KK=1
 1100       DO 1150 J=1,45
               DO 1140 JJ=1,3
                  IF(TRIP(J,JJ).EQ.EDGE .AND. JJ .EQ. 1) THEN
                     IF(II.EQ.1) THEN
                        PREV(KK,1)=TRIP(J,3)
                        PREV(KK,2)=TRIP(J,2)
                     ELSE
                        PREV(KK,1)=TRIP(J,2)
                        PREV(KK,2)=TRIP(J,3)
                     END IF
                     KK=KK+1
                  ELSE IF (TRIP(J,JJ).EQ.EDGE .AND. JJ.EQ.2) THEN
                     IF(II.EQ.1) THEN
                        PREV(KK,1)=TRIP(J,1)
                        PREV(KK,2)=TRIP(J,3)
                     ELSE
                        PREV(KK,1)=TRIP(J,3)
                        PREV(KK,2)=TRIP(J,1)
                     END IF
                     KK=KK+1
                  ELSE IF(TRIP(J,JJ).EQ.EDGE) THEN
                     IF(II.EQ.1) THEN
                        PREV(KK,1)=TRIP(J,2)
			PREV(KK,2)=TRIP(J,1)
                     ELSE
                        PREV(KK,1)=TRIP(J,1)
			PREV(KK,2)=TRIP(J,2)
                     END IF
                     KK=KK+1
                  END IF
 1140          CONTINUE
 1150       CONTINUE


C     List the vertices that will be removed

            REMLIST(REMNB+1)=90*(NB-1)+EDGESLOT
            REMLIST(REMNB+2)=90*(NB-1)+2*PREV(1,1)-MOD(II,2)
            REMLIST(REMNB+3)=90*(NB-1)+2*PREV(2,1)-MOD(II,2)
            REMLIST(REMNB+4)=90*(NB-1)+2*PREV(3,1)-MOD(II,2)
            REMNB=REMNB+4

C            WRITE(*,*) REMLIST


C     Then add the edges from this "knot" to the "knots" that are
C     already corrected.

C Find the triples (PREV,*,*) different from (PREV, I,*) for x_i 
C Find the triples (PREV,*,*) different from (I, PREV,*) for y_i 

 1399       DO 1412 R=1,3
               L=1
 1400          DO 1411 J=1,45      
 1402             DO 1410 K=1,3
                     IF(II.EQ.1) THEN
                        IF(K.EQ.1) THEN
                           IF(TRIP(J,K).EQ.PREV(R,1) .AND. 
     +                          TRIP(J,K+1).NE.EDGE) THEN
                              ADDEDGE(R,L,1)=TRIP(J,K+1)
                              ADDEDGE(R,L,2)=TRIP(J,K+2)
                              L=L+1
                           END IF
                        ELSE IF(K.EQ.2) THEN
                           IF(TRIP(J,K).EQ.PREV(R,1) .AND. 
     +                          TRIP(J,K+1).NE.EDGE) THEN
                              ADDEDGE(R,L,1)=TRIP(J,K+1)
                              ADDEDGE(R,L,2)=TRIP(J,1)
                              L=L+1
                           END IF
                        ELSE                              
                           IF(TRIP(J,K).EQ.PREV(R,1) .AND. 
     +                          TRIP(J,1).NE.EDGE) THEN
                              ADDEDGE(R,L,1)=TRIP(J,1)
                              ADDEDGE(R,L,2)=TRIP(J,2)
                              L=L+1
                           END IF
                        END IF
                     ELSE
                        IF(K.EQ.1) THEN
                           IF(TRIP(J,K).EQ.PREV(R,1) .AND. 
     +                          TRIP(J,3).NE.EDGE) THEN
                              ADDEDGE(R,L,1)=TRIP(J,3)
                              ADDEDGE(R,L,2)=TRIP(J,2)
                              L=L+1
                           END IF
                        ELSE IF(K.EQ.2) THEN
                           IF(TRIP(J,K).EQ.PREV(R,1) .AND. 
     +                          TRIP(J,K-1).NE.EDGE) THEN
                              ADDEDGE(R,L,1)=TRIP(J,1)
                              ADDEDGE(R,L,2)=TRIP(J,3)
                              L=L+1
                           END IF
                        ELSE                              
                           IF(TRIP(J,K).EQ.PREV(R,1) .AND. 
     +                          TRIP(J,K-1).NE.EDGE) THEN
                              ADDEDGE(R,L,1)=TRIP(J,2)
                              ADDEDGE(R,L,2)=TRIP(J,1)
                              L=L+1
                           END IF
                        END IF
                     END IF                     
 1410             CONTINUE
 1411          CONTINUE
 1412       CONTINUE
C Add the found edges
      

C            WRITE(*,*) (((ADDEDGE(R,J,M),M=1,2),J=1,2),R=1,3)

C TARKISTA etta FLOORien korjaukset -1, -2 on oikein

 1420       DO 1440 R=1,3
 1421          DO 1439 J=1,2
                  IF(II.EQ.1) THEN
                     ROW=90*(NB-1)+2*ADDEDGE(R,J,1)
C                     ROW=90*(NB-1)+2*ADDEDGE(R,J,1)-1
                     COL=90*(2*PREV(R,2)/3-1)
     +                    +2*ADDEDGE(R,J,2)-1
C                     WRITE(*,*) 'x', LABEL, R, J
C                     WRITE(*,*) ADDEDGE(R,J,1), ROW
C                     WRITE(*,*) ADDEDGE(R,J,2), COL
C                     WRITE(*,*) ADJTWO(ROW,COL)
                     ADJTWO(ROW,COL)=1
                     ADJTWO(COL,ROW)=1
                  ELSE
                     ROW=90*(NB-1)+2*ADDEDGE(R,J,1)-1
                     COL=90*(2*FLOOR(REAL(PREV(R,2))/3))
     +                    +2*ADDEDGE(R,J,2)
C                     WRITE(*,*) 'y', LABEL
C                     WRITE(*,*) ADDEDGE(R,J,1), ROW
C                     WRITE(*,*) ADDEDGE(R,J,2), COL
C                     WRITE(*,*) ADJTWO(ROW,COL)
C                     WRITE(*,*) ROW
C                     WRITE(*,*) COL
                     ADJTWO(ROW,COL)=1
                     ADJTWO(COL,ROW)=1
                  END IF
 1439          CONTINUE
 1440       CONTINUE

C      WRITE(*,*) 'Part of the 2-link:'
C      WRITE(*,'(90I1)') ((ADJTWO(I,J), J=1,90),I=1,90)


 1499    CONTINUE
 1500 CONTINUE


C     Collecting the non-zero rows to FINAL


      KK=1

 1600 DO 1800 ROW=1,30*90

 1610    DO 1700 COL=1,30*90

            IF(ADJTWO(ROW,COL) .EQ. 1) THEN
               OK=1
               DO 1650 J=1,30*4
                  IF (ROW .EQ. REMLIST(J)) THEN
                     OK=0
                  END IF
                  IF (COL .EQ. REMLIST(J)) THEN
                     OK=0
                  END IF
 1650          CONTINUE
               IF(OK.EQ.1) THEN
C                  WRITE(*,*) ROW, COL, KK
                  FDATA(KK,1)=ROW
                  FDATA(KK,2)=COL
                  KK=KK+1
               END IF
            END IF

 1700    CONTINUE

 1800 CONTINUE

C          WRITE(*,*) FDATA(1,1)
C      WRITE(*,'(2I5)') ((FDATA(I,J),J=1,2),I=1,45)

      KK=1
      DO 2000 ROW=1,2700
         IF(FDATA(ROW,1).LT. FDATA(ROW+1,1)) THEN
            LDATA(FDATA(ROW,1))=KK
            KK=KK+1
         ELSE IF(FDATA(ROW,1).NE.0 .AND. FDATA(ROW+1,1).EQ.0) THEN
             LDATA(FDATA(ROW,1))=KK
         END IF
 2000 CONTINUE

      DO 2011 ROW=1,2700
         IF(FDATA(ROW,1).NE.0) THEN
            FDATA(ROW,1)=LDATA(FDATA(ROW,1))
            FDATA(ROW,2)=LDATA(FDATA(ROW,2))
         END IF
 2011 CONTINUE


         
C      WRITE(*,*) KK

      IF (TRIPNB.LT.9) THEN
         WRITE(7,'(A,I1,AI4,A,I4,A)') 'P',TRIPNB+1,':= Graph< 780|{', 
     +        FDATA(1,1),
     +        ',',FDATA(1,2),'}'
      ELSEIF(TRIPNB.LT.99) THEN
         WRITE(7,'(A,I2,AI4,A,I4,A)') 'P',TRIPNB+1,':= Graph< 780|{', 
     +        FDATA(1,1),
     +        ',',FDATA(1,2),'}' 
      ELSEIF(TRIPNB.LT.999) THEN
         WRITE(7,'(A,I3,AI4,A,I4,A)') 'P',TRIPNB+1,':= Graph< 780|{', 
     +        FDATA(1,1),
     +        ',',FDATA(1,2),'}' 
      ELSE
         WRITE(7,'(A,I4,AI4,A,I4,A)') 'P',TRIPNB+1,':= Graph< 780|{', 
     +        FDATA(1,1),
     +        ',',FDATA(1,2),'}' 
      END IF
      ROW=2
 3000 IF(FDATA(ROW+1,1).NE.0) THEN
         WRITE(7,'(A,I4,A,I4,A)') ',{',FDATA(ROW,1),',',FDATA(ROW,2),'}'
         ROW=ROW+1
         GOTO 3000
      END IF
      WRITE(7,'(A,I4,A,I4,A)') ',{',FDATA(ROW,1),',',FDATA(ROW,2),'}>;'

      IF(TRIPNB .GT.1) THEN
         IF (TRIPNB.LT.9) THEN
            WRITE(7,'(A,I1,A)') 'if IsIsomorphic(P',TRIPNB+1,',P1) then'
            WRITE(7,'(A)') 'print 1; end if;' 
            WRITE(7,'(A,I1,A)') 'if IsIsomorphic(P',TRIPNB+1,',P2) then'
            WRITE(7,'(A)') 'print 2; end if;' 
         ELSEIF(TRIPNB.LT.99) THEN
            WRITE(7,'(A,I2,A)') 'if IsIsomorphic(P',TRIPNB+1,',P1) then'
            WRITE(7,'(A)') 'print 1; end if;' 
            WRITE(7,'(A,I2,A)') 'if IsIsomorphic(P',TRIPNB+1,',P2) then'
            WRITE(7,'(A)') 'print 2; end if;' 
         ELSEIF(TRIPNB.LT.999) THEN
            WRITE(7,'(A,I3,A)') 'if IsIsomorphic(P',TRIPNB+1,',P1) then'
            WRITE(7,'(A)') 'print 1; end if;' 
            WRITE(7,'(A,I3,A)') 'if IsIsomorphic(P',TRIPNB+1,',P2) then'
            WRITE(7,'(A)') 'print 2; end if;' 
         ELSE
            WRITE(7,'(A,I4,A)') 'if IsIsomorphic(P',TRIPNB+1,',P1) then'
            WRITE(7,'(A)') 'print 1; end if;' 
            WRITE(7,'(A,I4,A)') 'if IsIsomorphic(P',TRIPNB+1,',P2) then'
            WRITE(7,'(A)') 'print 2; end if;' 
         END IF
      END IF

      IF(TRIPNB.LT. MAX-1) GOTO 100

 9000 FORMAT (26(I1,A))
 9999 END
