      SUBROUTINE CTIP (A, N1, N2, MOVED, NWORK, NDIM)
C ----------------------------------------------------------
C TRANSPOSITION OF A RECTANGULAR MATRIX IN SITU.
C BY NORMAN BRENNER, MIT, 1/72.  CF. ALG. 380, CACM, 5/70.
C TRANSPOSITION OF THE N1 BY N2 MATRIX A AMOUNTS TO
C REPLACING THE ELEMENT AT VECTOR POSITION I (0-ORIGIN)
C WITH THE ELEMENT AT POSITION N1*I (MOD N1*N2-1).
C EACH SUBCYCLE OF THIS PERMUTATION IS COMPLETED IN ORDER.
C ----------------------------------------------------------
      COMPLEX A(1)
      COMPLEX ATEMP, BTEMP
      INTEGER MOVED(NWORK)
      INTEGER IFACT(8), IPOWER(8), NEXP(8), IEXP(8)
      IF (N1.LT.2 .OR. N2.LT.2) GO TO 200
      N12 = N1*N2
      N = N1
      M = N12 - 1
      IF (N1.NE.N2) GO TO 30
C SQUARE MATRICES ARE DONE SEPARATELY FOR SPEED
      I1MIN = 2
      DO 20 I1MAX=N,M,N
        I2 = I1MIN + N - 1
        DO 10 I1=I1MIN,I1MAX
          ATEMP = A(I1)
          A(I1) = A(I2)
          A(I2) = ATEMP
          I2 = I2 + N
   10   CONTINUE
        I1MIN = I1MIN + N + 1
   20 CONTINUE
      RETURN
C MODULUS M IS FACTORED INTO PRIME POWERS.  EIGHT FACTORS
C SUFFICE UP TO M = 2*3*5*7*11*13*17*19 = 9,767,520.
   30 NDIM = 0
      CALL INFCTR(M, IFACT, IPOWER, NEXP, NPOWER)
      DO 40 IP=1,NPOWER
        IEXP(IP) = 0
   40 CONTINUE
C GENERATE EVERY DIVISOR OF M LESS THAN M/2
      IDIV = 1
      MHALF = M/2
   50 IF (IDIV.GE.MHALF) RETURN
C THE NUMBER OF ELEMENTS WHOSE INDEX IS DIVISIBLE BY IDIV
C AND BY NO OTHER DIVISOR OF M IS THE EULER TOTIENT
C FUNCTION, PHI(M/IDIV).
      NCOUNT = M/IDIV
      DO 60 IP=1,NPOWER
        IF (IEXP(IP).EQ.NEXP(IP)) GO TO 60
        NCOUNT = (NCOUNT/IFACT(IP))*(IFACT(IP)-1)
   60 CONTINUE
      IF (NWORK.LE.0) GO TO 75
      DO 70 I=1,NWORK
        MOVED(I) = 0
   70 CONTINUE
   75 ISTART = IDIV
C THE STARTING POINT OF A SUBCYCLE IS DIVISIBLE ONLY BY IDIV
C AND MUST NOT APPEAR IN ANY OTHER SUBCYCLE.
   80 MMIST = M - ISTART
      IF (ISTART.EQ.IDIV) GO TO 120
      NDIM = MAX0(NDIM,ISTART)
      IF (ISTART.GT.NWORK) GO TO 90
      IF (MOVED(ISTART).NE.0) GO TO 160
   90 ISOID = ISTART/IDIV
      DO 100 IP=1,NPOWER
        IF (IEXP(IP).EQ.NEXP(IP)) GO TO 100
        IF (MOD(ISOID,IFACT(IP)).EQ.0) GO TO 160
  100 CONTINUE
      IF (ISTART.LE.NWORK) GO TO 120
      ITEST = ISTART
  110 ITEST = MOD(N*ITEST,M)
      IF (ITEST.LT.ISTART .OR. ITEST.GT.MMIST) GO TO 160
      IF (ITEST.GT.ISTART .AND. ITEST.LT.MMIST) GO TO 110
  120 ATEMP = A(ISTART+1)
      BTEMP = A(MMIST+1)
      IA1 = ISTART
  130 IA2 = MOD(N*IA1,M)
      MMIA1 = M - IA1
      MMIA2 = M - IA2
      IF (IA1.LE.NWORK) MOVED(IA1) = 1
      IF (MMIA1.LE.NWORK) MOVED(MMIA1) = 1
      NCOUNT = NCOUNT - 2
C MOVE TWO ELEMENTS, THE SECOND FROM THE NEGATIVE
C SUBCYCLE.  CHECK FIRST FOR SUBCYCLE CLOSURE.
      IF (IA2.EQ.ISTART) GO TO 140
      IF (MMIA2.EQ.ISTART) GO TO 150
      A(IA1+1) = A(IA2+1)
      A(MMIA1+1) = A(MMIA2+1)
      IA1 = IA2
      GO TO 130
  140 A(IA1+1) = ATEMP
      A(MMIA1+1) = BTEMP
      GO TO 160
  150 A(IA1+1) = BTEMP
      A(MMIA1+1) = ATEMP
  160 ISTART = ISTART + IDIV
      IF (NCOUNT.GT.0) GO TO 80
      DO 180 IP=1,NPOWER
        IF (IEXP(IP).EQ.NEXP(IP)) GO TO 170
        IEXP(IP) = IEXP(IP) + 1
        IDIV = IDIV*IFACT(IP)
        GO TO 50
  170   IEXP(IP) = 0
        IDIV = IDIV/IPOWER(IP)
  180 CONTINUE
      RETURN
  200 IF (N1.NE.N2) NDIM = 0
      RETURN
      END
