C     ALGORITHM 641 COLLECTED ALGORITHMS FROM ACM.
C     ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.12, NO. 2,
C     JUN., 1986, P. 149.
C     PROGRAM TEST
C****************************************************************
C                                                               *
C    TEST EXAMPLES                                              *
C    =============                                              *
C                                                               *
C****************************************************************
C
      INTEGER AB(540),P(144),V(30),CHG(12),XM(45),NCOL(15),NCOL1(15),
     .        NROW(30),CONS(3),MK(2565),X0(2025),CODE(45),DIM,R,TR
      REAL NORM(30)
      INTEGER A1(8,5),B1(8,2),R0,RL(30,12),RR(12,15),A2(30,15),B2(30,3),
     .        A3(7,6),B3(7,2)
C
C****************************************************************
C                                                               *
C    EXAMPLE 1 :                                                *
C    MATRIX USED BY RAO, SUBRAMANIAN AND KRISHNAMURTHY (1976)   *
C                                                               *
C****************************************************************
C
      DATA A1,B1/22,14,-1,-3,9,9,2,4,10,7,13,-2,8,1,-6,5,2,10,-1,13,1,
     .     -7,6,0,3,0,-11,-2,-2,5,5,-2,7,8,3,4,4,-1,1,2,12,7,-14,-1,1,8,
     .     8,-1,1,0,0,0,0,0,0,0/
      M = 8
      N = 5
      NRS = 2
      DIM = 5
      NNRS = N + NRS
C
C   **PRINT THE MATRICES A1 AND B1
C
      WRITE (6,9001)
      DO 10 I = 1,M
         WRITE (6,9011) I, (A1(I,J),J=1,N)
   10 CONTINUE
      WRITE (6,9021)
      DO 20 I = 1,M
         WRITE (6,9011) I, (B1(I,J),J=1,NRS)
   20 CONTINUE
C
C   **COMPUTE MINPRM, MAXPRM AND TR.
C
      CALL NBTERM(M,N,NRS,A1,B1,NORM,MAXPRM,MINPRM,TR)
      WRITE (6,9031) MINPRM,MAXPRM,TR
C
C   **SOLVE THE SYSTEM OF LINEAR EQUATIONS A1*X=B1.
C
      MODE = 3
      CALL DRIVER(M,N,NRS,NNRS,DIM,A1,B1,MINPRM,MAXPRM,TR,MODE,AB,P,V,
     .            CHG,XM,R,NCOL,NCOL1,NROW,CONS,MK,X0,CODE)
C
C****************************************************************
C                                                               *
C    EXAMPLE 2 :                                                *
C    RANDOMLY GENERATED RANK-DEFICIENT MATRIX                   *
C                                                               *
C****************************************************************
C
      M = 30
      N = 15
      NRS = 3
      DIM = 12
      NNRS = N + NRS
C
C   **GENERATE THE MATRICES RL, RR AND B2 RANDOMLY.
C   **(ELEMENTS FROM -617 TO 618)
C
      R0 = 12
      KTEMP = 561
      DO 40 I = 1,M
         DO 30 J = 1,R0
            KTEMP = MOD(KTEMP*2543,1237)
            RL(I,J) = KTEMP - 618
   30    CONTINUE
   40 CONTINUE
      DO 60 I = 1,R0
         DO 50 J = 1,N
            KTEMP = MOD(KTEMP*2543,1237)
            RR(I,J) = KTEMP - 618
   50    CONTINUE
   60 CONTINUE
      DO 80 I = 1,M
         DO 70 J = 1,NRS
            KTEMP = MOD(KTEMP*2543,1237)
            B2(I,J) = KTEMP - 618
   70    CONTINUE
   80 CONTINUE
C
C   **COMPUTE A2=RL*RR.
C
      DO 110 I = 1,M
         DO 100 J = 1,N
            KTEMP = 0
            DO 90 K = 1,R0
               KTEMP = KTEMP + RL(I,K)*RR(K,J)
   90       CONTINUE
            A2(I,J) = KTEMP
  100    CONTINUE
  110 CONTINUE
C
C   **PRINT THE MATRICES A2 AND B2.
C
      WRITE (6,9001)
      DO 120 I = 1,M
         WRITE (6,9011) I, (A2(I,J),J=1,N)
  120 CONTINUE
      WRITE (6,9021)
      DO 130 I = 1,M
         WRITE (6,9011) I, (B2(I,J),J=1,NRS)
  130 CONTINUE
C
C   **COMPUTE MINPRM, MAXPRM AND TR.
C
      CALL NBTERM(M,N,NRS,A2,B2,NORM,MAXPRM,MINPRM,TR)
      WRITE (6,9031) MINPRM,MAXPRM,TR
C
C   **SOLVE THE SYSTEM A2*X=B2.
C
      CALL DRIVER(M,N,NRS,NNRS,DIM,A2,B2,MINPRM,MAXPRM,TR,MODE,AB,P,V,
     .            CHG,XM,R,NCOL,NCOL1,NROW,CONS,MK,X0,CODE)
C
C****************************************************************
C                                                               *
C    EXAMPLE 3 :                                                *
C    ILL-CONDITIONED TEST MATRIX USED BY ZIELKE (1977)          *
C    (PARAMETER  A = 100000)                                    *
C                                                               *
C****************************************************************
C
      M = 7
      N = 6
      NRS = 2
      DIM = 6
      NNRS = N + NRS
C
C   **GENERATE THE MATRICES A3 AND B3.
C
      DO 150 K = 1,6
         IENTRY = 100007 - K
         DO 140 I = 1,K
            A3(I,K) = IENTRY
            A3(K,I) = IENTRY
  140    CONTINUE
  150 CONTINUE
      A3(5,5) = 100001
      A3(5,6) = 100000
      A3(6,5) = 100000
      A3(6,6) = 99999
      DO 160 J = 1,4
         A3(7,J) = 100000
  160 CONTINUE
      A3(7,5) = 99999
      A3(7,6) = 99998
      DO 170 I = 1,M
         B3(I,1) = 1
         B3(I,2) = I
  170 CONTINUE
C
C   **PRINT THE MATRICES A3 AND B3.
C
      WRITE (6,9001)
      DO 180 I = 1,M
         WRITE (6,9011) I, (A3(I,J),J=1,N)
  180 CONTINUE
      WRITE (6,9021)
      DO 190 I = 1,M
         WRITE (6,9011) I, (B3(I,J),J=1,NRS)
  190 CONTINUE
C
C   **COMPUTE MINPRM, MAXPRM AND TR.
C
      CALL NBTERM(M,N,NRS,A3,B3,NORM,MAXPRM,MINPRM,TR)
      WRITE (6,9031) MINPRM,MAXPRM,TR
C
C   **SOLVE THE SYSTEM OF LINEAR EQUATIONS A3*X=B3.
C
      CALL DRIVER(M,N,NRS,NNRS,DIM,A3,B3,MINPRM,MAXPRM,TR,MODE,AB,P,V,
     .            CHG,XM,R,NCOL,NCOL1,NROW,CONS,MK,X0,CODE)
C
      STOP

 9001 FORMAT (1H1,16X,6HINPUTS/17X,6 (1H*)//6X,
     .  24HMATRIX A OF COEFFICIENTS/)
 9011 FORMAT (5X,I3,1H),3X,9I12/ (12X,9I12))
 9021 FORMAT (/6X,28HMATRIX B OF RIGHT-HAND SIDES/)
 9031 FORMAT (/6X,9HMINPRM = ,I3,5X,9HMAXPRM = ,I3,5X,5HTR = ,I3/)
      END
C     SUBROUTINE NBTERM(M,N,NRS,A,B,NORM,MAXPRM,MINPRM,TR)
C
C****************************************************************
C                                                               *
C USING SEVERAL INEQUALITIES DESCRIBED IN THE COMPANION PAPER   *
C THIS SUBROUTINE COMPUTES THREE INTEGER NUMBERS MAXPRM, MINPRM *
C AND TR WHICH ARE USED BY EXSOLG TO TERMINATE THE EXACT        *
C SOLUTION OF THE SYSTEM OF LINEAR EQUATIONS AX=B.              *
C                                                               *
C  PRIME  IS A LINEAR ARRAY CONTAINING 100 DISTINCT PRIME       *
C (INPUT) INTEGERS IN ASCENDING ORDER. THE PRIMES ARE CHOSEN    *
C         AS LARGE AS POSSIBLE SUBJECT TO THE CONDITION THAT    *
C         FOR ALL I AND J PRIME(I)*PRIME(J) DOES NOT OVERFLOW   *
C         AN INTEGER WORD. THESE PRIMES ARE USED BY EXSOLG      *
C         AS MODULI IN THE EXACT COMPUTATION AND AS RADII FOR   *
C         THE REPRESENTATION OF THE INTEGER RESULTS IN          *
C         MIXED-RADIX FROM.                                     *
C      M  IS THE NUMBER OF EQUATIONS IN THE SYSTEM. (I.E., M IS *
C (INPUT) THE SIZE OF THE FIRST DIMENSIONS OF A AND B.)         *
C      N  IS THE NUMBER OF UNKNOWNS IN THE SYSTEM. (I.E., N IS  *
C (INPUT) THE SIZE OF THE SECOND DIMENSION OF A.)               *
C    NRS  IS THE NUMBER OF RIGHT-HAND SIDES IN THE SYSTEM.      *
C (INPUT) (I.E., NRS IS THE SIZE OF THE SECOND DIMENSION OF B.) *
C      A  IS AN INTEGER MATRIX OF DIMENSION M BY N WHICH        *
C (INPUT) CONTAINS THE MATRIX OF COEFFICIENTS.                  *
C      B  IS AN INTEGER MATRIX OF DIMENSION M BY NRS WHICH      *
C (INPUT) CONTAINS THE MATRIX OF RIGHT-HAND SIDES.              *
C   NORM  IS A REAL ARRAY OF DIMENSION MAX(M,N) USED TO STORE   *
C  (TEMP) THE EUCLIDEAN NORMS OF THE ROWS (IF M>=N) OR COLUMNS  *
C         (IF M<N) OF A                                         *
C MAXPRM  IS A NUMBER OF PRIMES SUFFICIENT FOR THE EXACT        *
C(OUTPUT) SOLUTION OF THE SYSTEM.                               *
C MINPRM  IS A NUMBER OF PRIMES SUFFICIENT FOR THE CORRECT      *
C(OUTPUT) RANK-DETERMINATION AND THE CORRECT DECISION ON        *
C         CONSISTENCY.                                          *
C     TR  IS THE MINIMUM NUMBER OF CONSECUTIVE ZERO COEFFICIENTS*
C(OUTPUT) IN THE MIXED-RADIX FORM OF THE INTEGER RESULTS OF     *
C         EXSOLG REQUIRED FOR THE TERMINATION BEFORE MAXPRM     *
C         PRIMES WERE USED.                                     *
C                                                               *
C****************************************************************
C
      SUBROUTINE NBTERM(M,N,NRS,A,B,NORM,MAXPRM,MINPRM,TR)
      INTEGER A(M,N),B(M,NRS),R,TR,PRIME(100)
      REAL NORM(1),MAX,MAX1,MAX2
      COMMON /PRIMEB/PRIME,IPRIME(100)
C
C   **COMPUTE  MAX = MAX(A(I,J)**2,B(I,J)**2) ,
C   **       FNORM = (FROBENIUS NORM OF A)**2
C   **AND THE SQUARES OF THE EUCLIDEAN NORMS OF THE COLUMNS OR ROWS
C   **OF A.
C
      MAX = 0.
      FNORM = 0.
      IF (M.LT.N) GO TO 10
      GO TO 40
C     THEN
   10 R = M
      K = N
      DO 30 J = 1,N
         S = 0.
         DO 20 I = 1,M
            X = FLOAT(A(I,J))**2
            FNORM = FNORM + X
            S = S + X
            IF (X.GT.MAX) MAX = X
   20    CONTINUE
         NORM(J) = S
   30 CONTINUE
      GO TO 70
C     ELSE
   40 R = N
      K = M
      DO 60 I = 1,M
         S = 0.
         DO 50 J = 1,N
            X = FLOAT(A(I,J))**2
            FNORM = FNORM + X
            S = S + X
            IF (X.GT.MAX) MAX = X
   50    CONTINUE
         NORM(I) = S
   60 CONTINUE
   70 RR = FLOAT(R)
      FNORM = AMAX1(FNORM,RR)
C
C   **COMPUTE  MAX1 = MAXIMUM SUM OF ABSOLUTE VALUES IN A COLUMN OF B,
C   **         MAX2 = MAXIMUM SUM OF SQUARES IN A COLUMN OF B.
C
      MAX1 = 1.
      MAX2 = 1.
      DO 90 J = 1,NRS
         S1 = 0.
         S2 = 0.
         DO 80 I = 1,M
            X = FLOAT(IABS(B(I,J)))
            S1 = S1 + X
            X = X**2
            IF (X.GT.MAX) MAX = X
            S2 = S2 + X
   80    CONTINUE
         IF (S1.GT.MAX1) MAX1 = S1
         IF (S2.GT.MAX2) MAX2 = S2
   90 CONTINUE
C   **COMPUTE THE BOUND B1 FOR THE DETERMINATION OF MAXPRM.
C
      ALOG2 = ALOG(2.)
      B1 = RR*ALOG(FNORM/RR)
      RR1 = AMAX1(1.,FLOAT(R-1))
      X = 0.5* ((2.*RR-1.)*ALOG(FNORM)+ALOG(MAX2)-RR*ALOG(RR)-RR1*
     .    ALOG(RR1))
      IF (X.GT.B1) B1 = X
      B1 = B1 + ALOG2
C
C   **COMPUTE THE BOUND B3 FOR THE DETERMINATION OF TR.
C
      B3 = ALOG(0.5*FLOAT(M* (N+1))*MAX)
C
C   **COMPUTE THE BOUND B2 FOR THE DETERMINATION OF MINPRM.
C
      B2 = 0.
      DO 120 I = 1,R
         MAX = 1.
         DO 110 J = 1,K
            IF (NORM(J).GT.MAX) GO TO 100
            GO TO 110
C           THEN
  100       MAX = NORM(J)
            L = J
  110    CONTINUE
         IF (MAX.LE.1.) GO TO 130
C        THEN
C     <-----EXIT
C        ELSE
         B2 = B2 + ALOG(MAX)
         NORM(L) = 0.
  120 CONTINUE
  130 B2 = 0.5*B2 + ALOG(MAX1) + ALOG2
C
C   **DETERMINE THE NUMBERS MAXPRM, MINPRM AND TR.
C
      MAXPRM = 0
      MINPRM = 1
      TR = 1
      SUMLOG = 0.
C     REPEAT
  140 MAXPRM = MAXPRM + 1
      SUMLOG = SUMLOG + ALOG(FLOAT(PRIME(MAXPRM)))
      IF (SUMLOG.LT.B2) MINPRM = MINPRM + 1
      IF (SUMLOG.LT.B3) TR = TR + 1
      IF (SUMLOG.GE.B1) RETURN
      GO TO 140
C        THEN
C     <-----RETURN
C     CONTINUE
      END
      SUBROUTINE DRIVER(M,N,NRS,NNRS,DIM,A,B,MINPRM,MAXPRM,TR,MODE,AB,P,
     .                  V,CHG,XM,R,NCOL,NCOL1,NROW,CONS,MK,X0,CODE)
C
C*****************************************************************
C                                                                *
C              CALL EXSOLG.                                      *
C              CALL FRADIX.                                      *
C              PRINT OUTPUTS IN FIXED-RADIX FORM.                *
C                                                                *
C*****************************************************************
C
      INTEGER DIM,TR,R
      INTEGER A(M,N),B(M,NRS),AB(M,NNRS),P(DIM,DIM),V(M),CHG(DIM),
     .        XM(N,NRS),NCOL(N),NCOL1(N),NROW(M),CONS(NRS),MK(MINPRM,1),
     .        X0(MAXPRM,N,NRS),CODE(MAXPRM),AP
      INTEGER A1(120),C1(120),BASE
      COMMON /ACTPR/PP,IP,P2,AP(100)
C
C   **SOLVE THE SYSTEM OF LINEAR EQUATIONS AX=B.
C
      CALL EXSOLG(M,N,NRS,NNRS,DIM,A,B,MINPRM,MAXPRM,TR,MODE,AB,P,V,CHG,
     .            XM,R,NCOL,NROW,CONS,MK,X0,CODE,NOCOEF,IER)
C
C   **PRINT OUTPUT PARAMETERS.
C
      WRITE (6,9001) IER,NOCOEF,R
      IF (R.NE.0) WRITE (6,9011) (NROW(I),I=1,R)
      IF (R.NE.0) WRITE (6,9021) (NCOL(I),I=1,R)
C
C   **CONVERT THE INTEGER RESULTS FROM MIXED-RADIX TO FIXED-RADIX FORM.
C   **PRINT THE INTEGER RESULTS IN FIXED-RADIX FORM.
C
      BASE = 10000
      LMAX = 120
      WRITE (6,9031)
      IF (R.EQ.0) GO TO 130
C     THEN
C      **NO INFORMATION ON THE NULL SPACE IS PRINTED.
      IF (R.EQ.N) GO TO 10
      GO TO 20
C     THEN
   10 WRITE (6,9041)
      GO TO 130
C     ELSE
   20 WRITE (6,9051)
      NMR = N - R
      K = R*NMR + 1
      DO 30 I = 1,MINPRM
         C1(I) = -MK(I,K)
   30 CONTINUE
      CALL FRADIX(C1,MINPRM,A1,LMAX,L0,BASE,IER)
      DO 60 I = 1,N
C         **DETERMINE II WITH NCOL(II)=I AND LET NCOL1(I)=II.
         DO 40 II = 1,N
            IF (NCOL(II).EQ.I) GO TO 50
C              THEN
C           <-----EXIT
   40    CONTINUE
   50    NCOL1(I) = II
   60 CONTINUE
      DO 120 J = 1,NMR
         WRITE (6,9061) J
         DO 110 I = 1,N
            II = NCOL1(I)
            IF (II.LE.R) GO TO 70
            GO TO 80
C              THEN
   70       K = (J-1)*R + II
            CALL FRADIX(MK(1,K),MINPRM,C1,LMAX,L,BASE,IER)
            WRITE (6,9071) (C1(K),K=1,L)
            GO TO 110
C              ELSE
   80       II = II - R
            IF (II.EQ.J) GO TO 90
            GO TO 100
C                 THEN
   90       WRITE (6,9071) (A1(K),K=1,L0)
            GO TO 110
C                 ELSE
  100       WRITE (6,9081)
  110    CONTINUE
  120 CONTINUE
C
  130 WRITE (6,9091)
      IF (MODE.EQ.1) WRITE (6,9101)
      IF (MODE.EQ.2) WRITE (6,9111)
      IF (MODE.EQ.3) WRITE (6,9121)
      WRITE (6,9131)
      DO 150 J = 1,NRS
         WRITE (6,9141) J
         IF (CONS(J).EQ.1) WRITE (6,9151)
         IF (CONS(J).EQ.0) WRITE (6,9161)
         WRITE (6,9171)
         DO 140 I = 1,N
            CALL FRADIX(X0(1,I,J),NOCOEF,C1,LMAX,L,BASE,IER)
            WRITE (6,9071) (C1(K),K=1,L)
  140    CONTINUE
  150 CONTINUE
      WRITE (6,9181)
      CALL FRADIX(CODE,NOCOEF,C1,LMAX,L,BASE,IER)
      WRITE (6,9071) (C1(K),K=1,L)
      RETURN

 9001 FORMAT (/17X,7HOUTPUTS/17X,7 (1H*)//6X,6HIER = ,I3,5X,9HNOCOEF = ,
     .  I3//6X,42HTHE MATRIX A OF COEFFICIENTS HAS RANK R = ,I3,1H.)
 9011 FORMAT (6X,44HTHE ROW SPACE OF A IS SPANNED BY THE LINEAR ,
     .  16HINDEPENDENT ROWS/8X,25 (I3,1H,))
 9021 FORMAT (6X,47HTHE COLUMN SPACE OF A IS SPANNED BY THE LINEAR ,
     .  19HINDEPENDENT COLUMNS/8X,25 (I3,1H,))
 9031 FORMAT (//6X,92 (1H=)/6X,36HEACH OF THE INTEGERS GIVEN BELOW IS ,
     .          53HIN FIXED-RADIX FORM. FROM LEFT TO RIGHT, THE NUMBERS
     .                    ,4HARE /6X,
     .                    43HCOEF(1),COEF(2),...,COEF(L), AND THE VALUE
     .                    ,17HOF THE INTEGER IS/10X,
     .                    20HCOEF(1)*BASE**(L-1)+,
     .                    31HCOEF(2)*BASE**(L-2)+...+COEF(L)/6X,
     .                    18HWITH BASE = 10000./6X,92 (1H=))
 9041 FORMAT (//6X,43HTHE NULL VECTOR IS THE ONLY ELEMENT OF THE ,
     .  16HNULL SPACE OF A.)
 9051 FORMAT (//6X,49HTHE NULL SPACE OF A IS SPANNED BY THE SUBSEQUENT ,
     .  34HLINEAR INDEPENDENT INTEGER VECTORS)
 9061 FORMAT (/10X,11HVECTOR NO. ,I3/)
 9071 FORMAT (14X,17I6/ (20X,16I6))
 9081 FORMAT (19X,1H0)
 9091 FORMAT (//6X,38HFOR EACH RIGHT-HAND SIDE OF THE SYSTEM)
 9101 FORMAT (6X,46HA PARTICULAR SOLUTION OF R LINEAR INDEPENDENT ,
     .  23HEQUATIONS OF THE SYSTEM)
 9111 FORMAT (6X,52HA PARTICULAR SOLUTION OR A PARTICULAR LEAST SQUARES
     .  ,8HSOLUTION)
 9121 FORMAT (6X,44HTHE SOLUTION OR LEAST SQUARES SOLUTION WITH ,
     .  22HMINIMUM EUCLIDEAN NORM)
 9131 FORMAT (6X,13HWAS COMPUTED.)
 9141 FORMAT (/10X,20HRIGHT-HAND SIDE NO. ,I3)
 9151 FORMAT (/12X,25HTHE SYSTEM IS CONSISTENT.)
 9161 FORMAT (/12X,27HTHE SYSTEM IS INCONSISTENT.)
 9171 FORMAT (12X,42HNUMERATORS OF THE ELEMENTS OF THE SOLUTION/)
 9181 FORMAT (/10X,36HCOMMON DENOMINATOR FOR ALL SOLUTIONS/)
      END
C     SUBROUTINE EXSOLG(M,N,NRS,NNRS,DIM,A,B,MINPRM,MAXPRM,TR,MODE,
C                       AB,P,V,CHG,XM,
C                       R,NCOL,NROW,CONS,MK,X0,CODE,NOCOEF,IER)
C
C****************************************************************
C                                                               *
C THIS SUBROUTINE SOLVES EXACTLY THE GENERAL SYSTEM OF LINEAR   *
C EQUATIONS AX=B, WHERE A(M,N) AND B(M,NRS) ARE MATRICES WITH   *
C INTEGER ENTRIES. IT COMPUTES THE RANK R, A SET OF R LINEAR    *
C INDEPENDENT ROWS AND COLUMNS AND THE MATRIX K FROM THE        *
C HERMITIAN STANDARD FORM OF A. IT DECIDES ON THE CONSISTENCY   *
C OF THE EQUATIONS FOR EACH RIGHT-HAND SIDE AND COMPUTES        *
C A SOLUTION CORRESPONDING TO THE VALUE OF MODE.                *
C THE COMMON DENOMINATOR OF THE ELEMENTS OF K AND THE           *
C CORRESPONDING NUMERATORS ARE STORED IN THE ARRAY MK. THE      *
C COMMON DENOMINATOR OF THE SOLUTION IS STORED IN CODE AND THE  *
C CORRESPONDING NUMERATORS ARE STORED IN X0. ALL THESE INTEGER  *
C VALUES ARE IN MIXED-RADIX FORM.                               *
C IF THE SOLUTION X IS REQUIRED, THE USER NEED ONLY COMPUTE     *
C X=X0/CODE.                                                    *
C                                                               *
C  PRIME  IS A LINEAR ARRAY CONTAINING 100 DISTINCT PRIME       *
C (INPUT) INTEGERS IN ASCENDING ORDER. THE PRIMES ARE CHOSEN    *
C         AS LARGE AS POSSIBLE SUBJECT TO THE CONDITION THAT    *
C         FOR ALL I AND J PRIME(I)*PRIME(J) DOES NOT OVERFLOW   *
C         AN INTEGER WORD. THESE PRIMES ARE USED AS RADII FOR   *
C         THE REPRESENTATION OF THE INTEGER RESULTS IN          *
C         MIXED-RADIX FORM.                                     *
C IPRIME  IS A LINEAR ARRAY OF INTEGERS SUCH THAT FOR EACH K,   *
C (INPUT) IPRIME(K)* PRIME(1)*PRIME(2)*...*PRIME(K-1) = 1       *
C                                          MODULO PRIME(K).     *
C     AP  IS A LINEAR INTEGER ARRAY OF DIMENSION 100 WHICH      *
C(OUTPUT) CONTAINS THE ACTUALLY USED PRIMES. THIS ARRAY IS      *
C         NECESSARY BECAUSE IN RARE CASES SOME PRIMES FROM THE  *
C         ARRAY PRIME ARE NOT SUITABLE FOR THE SPECIAL SYSTEM   *
C         OF EQUATIONS AND MUST BE DISCARDED.                   *
C      M  IS THE NUMBER OF EQUATIONS IN THE SYSTEM. (I.E.,      *
C (INPUT) M IS THE SIZE OF THE FIRST DIMENSIONS OF A AND B.)    *
C      N  IS THE NUMBER OF UNKNOWNS IN THE SYSTEM. (I.E., N IS  *
C (INPUT) THE SIZE OF THE SECOND DIMENSION OF A.)               *
C    NRS  IS THE NUMBER OF RIGHT-HAND SIDES FOR WHICH THE SYSTEM*
C (INPUT) IS TO BE SOLVED. (I.E., NRS IS THE SIZE OF THE SECOND *
C         DIMENSION OF B.)                                      *
C   NNRS  IS EQUAL TO N+NRS.                                    *
C (INPUT)                                                       *
C    DIM  IS THE SIZE OF THE FIRST AND SECOND DIMENSION OF THE  *
C (INPUT) MATRIX P. DIM MUST BE EQUAL TO OR GREATER THAN THE    *
C         RANK R OF A, MIN(M,N) IS ALWAYS SUFFICIENT.           *
C      A  IS AN INTEGER ARRAY OF DIMENSION M BY N WHICH CONTAINS*
C (INPUT) THE MATRIX OF COEFFICIENTS OF THE SYSTEM.             *
C      B  IS AN INTEGER ARRAY OF DIMENSION M BY NRS WHICH       *
C (INPUT) CONTAINS THE MATRIX OF RIGHT-HAND SIDES OF THE SYSTEM.*
C MINPRM  IS THE MINIMUM NUMBER OF PRIMES TO BE USED.           *
C (INPUT)                                                       *
C MAXPRM  IS THE MAXIMUM NUMBER OF PRIMES TO BE USED. IT MUST BE*
C (INPUT) LESS THEN 101.                                        *
C     TR  IS THE NUMBER OF CONSECUTIVE ZERO COEFFICIENTS IN THE *
C (INPUT) MIXED-RADIX REPRESENTATION OF X0 AND CODE REQUIRED TO *
C         GUARANTEE THAT X=X0/CODE IS THE (GENERALIZED) SOLUTION*
C         OF AX=B.                                              *
C         SUITABLE VALUES FOR MINPRM, MAXPRM AND TR CAN BE      *
C         OBTAINED USING SUBROUTINE NBTERM.                     *
C   MODE  IS A CONTROL PARAMETER WHICH DESCRIBES THE SOLUTION   *
C (INPUT) TO BE COMPUTED:                                       *
C         IF MODE = 1, A PARTICULAR SOLUTION OF R LINEAR        *
C            INDEPENDENT EQUATIONS IS COMPUTED. IT IS A         *
C            PARTICULAR SOLUTION OF THE COMPLETE SYSTEM IF THE  *
C            SYSTEM IS CONSISTENT.                              *
C         IF MODE = 2, A PARTICULAR SOLUTION (CONSISTENT CASE)  *
C            OR A PARTICULAR LEAST SQUARES SOLUTION             *
C            (INCONSISTENT CASE) IS COMPUTED.                   *
C         IF MODE = 3, THE SOLUTION OR LEAST SQUARES SOLUTION   *
C            WITH MINIMUM EUCLIDEAN NORM IS COMPUTED.           *
C     AB  IS AN INTEGER ARRAY OF DIMENSION M BY N+NRS USED FOR  *
C  (TEMP) TEMPORARY STORAGE OF THE AUGMENTED MATRIX (A,B) MODULO*
C         THE VARIOUS PRIMES PRIME(K).                          *
C      P  IS AN INTEGER ARRAY OF DIMENSION DIM BY DIM USED FOR  *
C  (TEMP) TEMPORARY STORAGE OF CERTAIN MATRIX PRODUCTS.         *
C      V  IS A LINEAR INTEGER ARRAY OF DIMENSION M USED FOR     *
C  (TEMP) TEMPORARY STORAGE OF COLUMNS.                         *
C    CHG  IS A LINEAR INTEGER ARRAY OF DIMENSION DIM USED FOR   *
C  (TEMP) DESCRIBING PERMUTATIONS OF ROWS AND COLUMNS BY THE    *
C         SUBROUTINE SYMSOL.                                    *
C     XM  IS AN INTEGER ARRAY OF DIMENSION N BY NRS USED FOR    *
C  (TEMP) TEMPORARY STORAGE OF THE SOLUTIONS MODULO THE VARIOUS *
C         PRIMES PRIME(K).                                      *
C         THE ARRAYS AB, P, V, CHG AND XM ARE INCLUDED IN THE   *
C         ARGUMENT LIST ONLY TO PERMIT THEIR DIMENSIONS TO BE   *
C         VARIABLE.                                             *
C      R  IS THE RANK OF THE MATRIX A OF COEFFICIENTS.          *
C(OUTPUT)                                                       *
C   NCOL  IS A LINEAR INTEGER ARRAY OF DIMENSION N WHICH        *
C(OUTPUT) CONTAINS THE INDICES OF THE COLUMNS OF A IN THE ORDER *
C         AS THEY ARE USED IN THE TRANSFORMATION. THE FIRST R   *
C         ELEMENTS OF NCOL GIVE THE INDICES OF R LINEAR         *
C         INDEPENDENT COLUMNS OF A.                             *
C   NROW  IS A LINEAR INTEGER ARRAY OF DIMENSION M WHICH        *
C(OUTPUT) CONTAINS THE INDICES OF THE ROWS OF A IN THE ORDER    *
C         AS THEY ARE USED IN THE TRANSFORMATION. THE FIRST R   *
C         ELEMENTS OF NROW GIVE THE INDICES OF R LINEAR         *
C         INDEPENDENT ROWS OF A.                                *
C   CONS  IS A LINEAR INTEGER ARRAY OF DIMENSION NRS.           *
C(OUTPUT) CONS(I) = 1 IF THE SYSTEM WITH THE I-TH RIGHT-HAND    *
C                     SIDE IS CONSISTENT.                       *
C         CONS(I) = 0 IF THE SYSTEM WITH THE I-TH RIGHT-HAND    *
C                     SIDE IS INCONSISTENT.                     *
C     MK  IS AN INTEGER ARRAY OF TWO DIMENSIONS. THE SIZE OF THE*
C(OUTPUT) FIRST DIMENSION IS MINPRM. THE SIZE OF THE SECOND     *
C         DIMENSION IS AT LEAST R*(N-R)+1, INT(N**2/4)+1 IS     *
C         ALWAYS SUFFICIENT. IF R IS NOT EQUAL N, MK CONTAINS   *
C         THE COEFFICIENTS OF THE MIXED-RADIX REPRESENTATIONS   *
C         OF THE NUMERATORS (ELEMENTS 1 TO R*(N-R), COLUMNWISE) *
C         AND OF THE COMMON DENOMINATOR (ELEMENT R*(N-R)+1) OF  *
C         THE ELEMENTS OF THE MATRIX K FROM THE HERMITIAN       *
C         STANDARD FORM:                                        *
C         NUMERATOR OF K(I,J)                                   *
C         =   MK(1,R*(J-1)+I)                                   *
C           + MK(2,R*(J-1)+I)*AP(1)                             *
C           .                                                   *
C           .                                                   *
C           .                                                   *
C           + MK(MINPRM,R*(J-1)+I)*AP(1)*AP(2)*...*AP(MINPRM-1) *
C         DENOMINATOR SIMILIAR.                                 *
C     X0  IS AN INTEGER ARRAY OF DIMENSION MAXPRM BY N BY NRS   *
C(OUTPUT) WHICH CONTAINS THE COEFFICIENTS OF THE MIXED-RADIX    *
C         REPRESENTATIONS OF THE NUMERATORS OF THE ELEMENTS OF  *
C         THE SOLUTION MATRIX:                                  *
C         X0( ,I,J)                                             *
C         =   X0(1,I,J)                                         *
C           + X0(2,I,J)*AP(1)                                   *
C           .                                                   *
C           .                                                   *
C           .                                                   *
C           + X0(NOCOEF,I,J)*AP(1)*AP(2)*...*AP(NOCOEF-1)       *
C   CODE  IS A LINEAR INTEGER ARRAY OF DIMENSION MAXPRM WHICH   *
C(OUTPUT) CONTAINS THE COEFFICIENTS OF THE MIXED-RADIX          *
C         REPRESENTATION OF THE COMMON DENOMINATOR OF THE       *
C         ELEMENTS OF THE SOLUTION MATRIX:                      *
C         CODE                                                  *
C         =   CODE(1)                                           *
C           + CODE(2)*AP(1)                                     *
C           .                                                   *
C           .                                                   *
C           .                                                   *
C           + CODE(NOCOEF)*AP(1)*AP(2)*...*AP(NOCOEF-1)         *
C NOCOEF  IS THE NUMBER OF PRIMES ACTUALLY USED TO REPRESENT    *
C(OUTPUT) X0( ,I,J) AND CODE IN MIXED-RADIX FORM.               *
C         NOCOEF <= MAXPRM.                                     *
C    IER  IS AN ERROR CODE WHICH IS                             *
C(OUTPUT) 0  IF THE SUBROUTINE HAS TERMINATED BY THE RECURSIVE  *
C            TEST.                                              *
C         1  IF THE SUBROUTINE HAS TERMINATED AFTER MAXPRM      *
C            OR MINPRM (IN THE CASE MODE=1) EXECUTIONS OF THE   *
C            BASIC ALGORITHM.                                   *
C         2  IF THE 100 PRIMES DID NOT SUFFICE.                 *
C         3  IF THE INPUT PARAMETERS ARE INCORRECT.             *
C         IF THE INPUT PARAMETERS MINPRM, MAXPRM AND TR WERE    *
C         SUITABLE CHOSEN (E.G. DETERMINED BY THE SUBROUTINE    *
C         NBTERM) AND IER <= 1, THE EXACT SOLUTION HAS BEEN     *
C         COMPUTED CORRECTLY.                                   *
C                                                               *
C****************************************************************
C
      SUBROUTINE EXSOLG(M,N,NRS,NNRS,DIM,A,B,MINPRM,MAXPRM,TR,MODE,AB,P,
     .                  V,CHG,XM,R,NCOL,NROW,CONS,MK,X0,CODE,NOCOEF,IER)
      INTEGER DIM,TR,R
      INTEGER A(M,N),B(M,NRS),AB(M,NNRS),P(DIM,DIM),V(M),CHG(DIM),
     .        XM(N,NRS),NCOL(N),NROW(M),CONS(NRS),MK(MINPRM,1),
     .        X0(MAXPRM,N,NRS),CODE(MAXPRM),AP,C,DIS,D0,D1,D2,PIV,PP,
     .        PRIME,P2,RP1,S,S1,SCAPRO
      LOGICAL TEST
      COMMON /PRIMEB/PRIME(100),IPRIME(100)
      COMMON /ACTPR/PP,IP,P2,AP(100)
      IER = 3
      IF (MAXPRM.GT.100 .OR. M.LT.1 .OR. N.LT.1 .OR. NRS.LT.1 .OR.
     .    NNRS.NE.N+NRS) RETURN
      IER = 2
      NP1 = N + 1
C
C   **DIS IS THE NUMBER OF DISCARDED PRIMES.
C
      DIS = 0
C
   10 R = 0
      DO 20 I = 1,NRS
         CONS(I) = 1
   20 CONTINUE
      DO 30 I = 1,N
         NCOL(I) = I
   30 CONTINUE
      DO 40 I = 1,M
         NROW(I) = I
   40 CONTINUE
C
C   **NZEROK IS THE NUMBER OF CONSECUTIVE ZERO COEFFICIENTS IN THE
C   **MIXED-RADIX REPRESENTATIONS OF THE DENOMINATOR AND THE NUMERATORS
C   **OF THE MATRIX K.
C
      NZEROK = 0
C
C   **NOZERO IS THE NUMBER OF CONSECUTIVE ZERO COEFFICIENTS IN THE
C   **MIXED-RADIX REPRESENTATIONS OF THE DENOMINATOR AND THE NUMERATORS
C   **OF THE COMPUTED SOLUTION.
C
      NOZERO = 0
C
      DO 1010 ITER = 1,MAXPRM
   50    ITDIS = ITER + DIS
         IF (ITDIS.GT.100) RETURN
         PP = PRIME(ITDIS)
C
C     ***THE SYSTEM AX=B IS SOLVED MODULO PP.
C
         P2 = (PP+1)/2
         ITER1 = ITER - 1
         ITER2 = ITER - 2
         D0 = 1
         D1 = 1
         D2 = 1
C
C      **AB IS THE AUGMENTED MATRIX (A,B) MODULO PP. THE ORDER OF THE
C      **COLUMNS OF A AND OF THE ROWS OF A AND B IS DETERMINED BY A
C      **PIVOT SEARCH FOR ITER=1.
C
         DO 80 I = 1,M
            II = NROW(I)
            DO 60 J = 1,N
               JJ = NCOL(J)
               AB(I,J) = MOD(A(II,JJ),PP)
   60       CONTINUE
            DO 70 J = 1,NRS
               JN = J + N
               AB(I,JN) = MOD(B(II,J),PP)
   70       CONTINUE
   80    CONTINUE
         S = 1
C
C        REPEAT ... TRANSFORMATION WITH THE PIVOT COLUMN S
C         **THE ELEMENTS OF THE TRANSFORMATION MATRIX ARE STORED
C         **IN THE ARRAY AB.
C
   90    IF (S.LE.R) GO TO 100
         GO TO 210
C           THEN
C            **TAKE THE SAME PIVOTS AS IN THE PREVIOUS EXECUTIONS.
C
  100    PIV = AB(S,S)
         IF (PIV.EQ.0) GO TO 110
         GO TO 120
C              THEN
C               **DROP THIS PRIME, GO BACK AND TAKE THE NEXT PRIME.
  110    DIS = DIS + 1
         GO TO 50

  120    D0 = MOD(D0*PIV,PP)
         PIV = INV(PIV,PP)
         K = 1
         IF (R.EQ.M .OR. MODE.EQ.1) K = S + 1
         DO 130 J = K,NNRS
            IF (J.NE.S) AB(S,J) = MOD(AB(S,J)*PIV,PP)
  130    CONTINUE
         AB(S,S) = PIV
         IF (S.GT.1) GO TO 140
         GO TO 170
C              THEN
  140    S1 = S - 1
         DO 160 I = 1,S1
            C = AB(I,S)
            DO 150 J = K,NNRS
               IF (J.NE.S) AB(I,J) = MOD(AB(I,J)-C*AB(S,J),PP)
  150       CONTINUE
            AB(I,S) = MOD(-C*PIV,PP)
  160    CONTINUE
  170    IF (S.EQ.M) GO TO 320
C              THEN
C        <--------EXIT
  180    S1 = S + 1
         DO 200 I = S1,M
            C = AB(I,S)
            DO 190 J = K,NNRS
               IF (J.NE.S) AB(I,J) = MOD(AB(I,J)-C*AB(S,J),PP)
  190       CONTINUE
            AB(I,S) = MOD(-C*PIV,PP)
  200    CONTINUE
         S = S1
         IF (S.GT.N) GO TO 320
C              THEN
C        <--------EXIT
         GO TO 90
C           ELSE
  210    IF (ITER.GT.MINPRM) GO TO 320
C              THEN
C        <--------EXIT
C
C            **LOOK FOR NONZERO PIVOTS.
C
         DO 310 I = S,M
            DO 300 J = S,N
               IF (AB(I,J).NE.0) GO TO 220
               GO TO 300
C                    THEN
C                     **MOVE AB(I,J) TO THE POSITION (S,S).
C
  220          IF (J.NE.S) GO TO 230
               GO TO 250
C                       THEN
  230          DO 240 II = 1,M
                  C = AB(II,S)
                  AB(II,S) = AB(II,J)
                  AB(II,J) = C
  240          CONTINUE
               C = NCOL(S)
               NCOL(S) = NCOL(J)
               NCOL(J) = C
  250          IF (I.NE.S) GO TO 260
               GO TO 280
C                       THEN
  260          DO 270 II = 1,NNRS
                  C = AB(S,II)
                  AB(S,II) = AB(I,II)
                  AB(I,II) = C
  270          CONTINUE
               C = NROW(S)
               NROW(S) = NROW(I)
               NROW(I) = C
  280          R = S
               IF (ITER.GT.1) GO TO 290
               GO TO 100
C                       THEN
C                        **DROP ALL PRIMES BEFORE AND GO BACK.
  290          DIS = DIS + ITER - 1
               GO TO 10

  300       CONTINUE
  310    CONTINUE
C        CONTINUE
C
C     ***TRANSFORMATION MODULO PP WAS SUCCESSFUL,
C     ***NOW COMPUTE THE SOLUTION.
C
  320    RP1 = R + 1
         AP(ITER) = PP
C
C      **COMPUTE IP WITH IP*AP(1)*...*AP(ITER-1)=1 MOD PP.
C
         IF (DIS.EQ.0 .OR. ITER.EQ.1) GO TO 330
         GO TO 340
C        THEN
  330    IP = IPRIME(ITER)
         GO TO 360
C        ELSE
  340    K = 1
         DO 350 I = 1,ITER1
            K = MOD(K*AP(I),PP)
  350    CONTINUE
         IP = INV(K,PP)
C
  360    IF (R.LT.M) GO TO 370
         GO TO 420
C        THEN
C         **CHECK FOR CONSISTENCY.
C
  370    DO 410 J = 1,NRS
            IF (CONS(J).EQ.1) GO TO 380
            GO TO 410
C              THEN
  380       JN = J + N
            DO 400 I = RP1,M
               IF (AB(I,JN).NE.0) GO TO 390
               GO TO 400
C                    THEN
  390          CONS(J) = 0
               GO TO 410

  400       CONTINUE
  410    CONTINUE
C
  420    IF (R.EQ.0) GO TO 890
C        THEN
C         **SET ALL ELEMENTS OF THE SOLUTION EQUAL TO ZERO
C         **AND TAKE THE NEXT PRIME.
C
         IF (ITER.LE.MINPRM .AND. R.LT.N) GO TO 430
         GO TO 480
C        THEN
C         **COMPUTE THE NEXT COEFFICIENTS OF THE MIXED-RADIX
C         **REPRESENTATIONS OF THE NUMERATORS AND THE DENOMINATOR
C         **OF THE MATRIK K.
C
  430    K = 0
         TEST = .TRUE.
         DO 450 J = RP1,N
            DO 440 I = 1,R
               K = K + 1
               CALL MIXRAD(MK(1,K),ITER,AB(I,J)*D0,TEST)
  440       CONTINUE
  450    CONTINUE
         K = K + 1
         CALL MIXRAD(MK(1,K),ITER,D0,TEST)
C
C         **COUNT THE CONSECUTIVE ZERO COEFFICIENTS.
C
         IF (TEST) GO TO 460
         GO TO 470
C           THEN
  460    NZEROK = NZEROK + 1
         GO TO 480
C           ELSE
  470    NZEROK = 0
C
  480    IF (MODE.EQ.1 .OR. R.EQ.M) GO TO 490
         GO TO 520
C        THEN
C         **STORE THE TRANSFORMED RIGHT-HAND SIDES IN XM.
C
  490    DO 510 I = 1,R
            DO 500 J = 1,NRS
               JN = J + N
               XM(I,J) = AB(I,JN)
  500       CONTINUE
  510    CONTINUE
         GO TO 710
C        ELSE
C         **COMPUTE THE PROJECTED AND TRANSFORMED RIGHT-HAND SIDES
C         **AND STORE THEM IN XM.
C
  520    MR = M - R
         IF (2*R.GE.M) GO TO 530
         GO TO 620
C           THEN
C            **USE FORMULA (A1).
C
  530    DO 550 I = 1,MR
            IR = I + R
            DO 540 J = I,MR
               JR = J + R
               P(I,J) = SCAPRO(AB(IR,1),M,AB(JR,1),M,R,PP)
  540       CONTINUE
            P(I,I) = MOD(P(I,I)+1,PP)
  550    CONTINUE
         CALL SYMSOL(P,DIM,MR,AB,M,RP1,NP1,NNRS,V,CHG,D1,PP,IE)
         IF (IE.NE.0) GO TO 560
         GO TO 570
C              THEN
C               **DROP THIS PRIME, GO BACK AND TAKE THE NEXT PRIME.
  560    DIS = DIS + 1
         GO TO 50

  570    DO 590 I = 1,R
            DO 580 J = 1,NRS
               JN = J + N
               XM(I,J) = SCAPRO(AB(RP1,I),I,AB(RP1,JN),1,MR,PP)
  580       CONTINUE
  590    CONTINUE
         DO 610 J = 1,NRS
            CALL VSTORE(XM(1,J),1,R,V)
            JN = J + N
            DO 600 I = 1,R
               C = SCAPRO(AB(I,1),M,V,1,R,PP)
               XM(I,J) = MOD(AB(I,JN)-C,PP)
  600       CONTINUE
  610    CONTINUE
         GO TO 710
C           ELSE
C            **USE FORMULA (A2).
C
  620    DO 640 I = 1,R
            DO 630 J = I,R
               P(I,J) = SCAPRO(AB(RP1,I),1,AB(RP1,J),1,MR,PP)
  630       CONTINUE
            P(I,I) = MOD(P(I,I)+1,PP)
  640    CONTINUE
         DO 660 J = NP1,NNRS
            CALL VSTORE(AB(RP1,J),1,MR,V)
            DO 650 I = 1,R
               IR = I + R
               AB(IR,J) = SCAPRO(AB(RP1,I),1,V,1,MR,PP)
  650       CONTINUE
  660    CONTINUE
         CALL SYMSOL(P,DIM,R,AB,M,RP1,NP1,NNRS,V,CHG,D1,PP,IE)
         IF (IE.NE.0) GO TO 670
         GO TO 680
C              THEN
C               **DROP THIS PRIME, GO BACK AND TAKE THE NEXT PRIME.
  670    DIS = DIS + 1
         GO TO 50

  680    DO 700 I = 1,R
            DO 690 J = 1,NRS
               JN = J + N
               C = SCAPRO(AB(I,1),M,AB(RP1,JN),1,R,PP)
               XM(I,J) = MOD(AB(I,JN)-C,PP)
  690       CONTINUE
  700    CONTINUE
C
  710    IF (MODE.EQ.3 .AND. R.LT.N) GO TO 720
         GO TO 890
C        THEN
C         **COMPUTE THE MINIMUM NORM SOLUTION
C
  720    NR = N - R
         IF (2*R.LT.N) GO TO 730
         GO TO 800
C           THEN
C            **USE FORMULA (B1).
C
  730    DO 750 I = 1,R
            DO 740 J = I,R
               AB(I,J) = SCAPRO(AB(I,RP1),M,AB(J,RP1),M,NR,PP)
  740       CONTINUE
            AB(I,I) = MOD(AB(I,I)+1,PP)
  750    CONTINUE
         CALL SYMSOL(AB,M,R,XM,N,1,1,NRS,V,CHG,D2,PP,IE)
         IF (IE.NE.0) GO TO 760
         GO TO 770
C              THEN
C               **DROP THIS PRIME, GO BACK AND TAKE THE NEXT PRIME
  760    DIS = DIS + 1
         GO TO 50

  770    DO 790 I = RP1,N
            DO 780 J = 1,NRS
               XM(I,J) = SCAPRO(AB(1,I),1,XM(1,J),1,R,PP)
  780       CONTINUE
  790    CONTINUE
         GO TO 890
C           ELSE
C            **USE FORMULA (B2).
C
  800    DO 820 I = 1,NR
            IR = I + R
            DO 810 J = I,NR
               JR = J + R
               AB(I,J) = SCAPRO(AB(1,IR),1,AB(1,JR),1,R,PP)
  810       CONTINUE
            AB(I,I) = MOD(AB(I,I)+1,PP)
  820    CONTINUE
         DO 840 I = RP1,N
            DO 830 J = 1,NRS
               XM(I,J) = SCAPRO(AB(1,I),1,XM(1,J),1,R,PP)
  830       CONTINUE
  840    CONTINUE
         CALL SYMSOL(AB,M,NR,XM,N,RP1,1,NRS,V,CHG,D2,PP,IE)
         IF (IE.NE.0) GO TO 850
         GO TO 860
C              THEN
C               **DROP THIS PRIME, GO BACK AND TAKE THE NEXT PRIME
  850    DIS = DIS + 1
         GO TO 50

  860    DO 880 I = 1,R
            DO 870 J = 1,NRS
               C = SCAPRO(AB(I,RP1),M,XM(RP1,J),1,NR,PP)
               XM(I,J) = MOD(XM(I,J)-C,PP)
  870       CONTINUE
  880    CONTINUE
C
  890    IF (MODE.GE.2 .AND. (R.NE.M.OR.R.NE.N)) D0 = MOD(D0*D0,PP)
         D0 = MOD(D0*MOD(D1*D2,PP),PP)
C
C      **COMPUTE THE NEXT COEFFICIENTS IN THE MIXED-RADIX
C      **REPRESENTATIONS OF THE NUMERATORS AND THE DENOMINATOR
C      **OF THE SOLUTION
C
         TEST = .TRUE.
         DO 940 I = 1,N
            II = NCOL(I)
            IF (I.GT.R .AND. MODE.LE.2 .OR. R.EQ.0) GO TO 900
            GO TO 920
C           THEN
  900       DO 910 J = 1,NRS
               X0(ITER,II,J) = 0
  910       CONTINUE
            GO TO 940
C           ELSE
  920       DO 930 J = 1,NRS
               CALL MIXRAD(X0(1,II,J),ITER,D0*XM(I,J),TEST)
  930       CONTINUE
  940    CONTINUE
         CALL MIXRAD(CODE,ITER,D0,TEST)
C
C      **COUNT THE CONSECUTIVE ZERO COEFFICIENTS
C
         IF (TEST) GO TO 950
         GO TO 960
C        THEN
  950    NOZERO = NOZERO + 1
         GO TO 970
C        ELSE
  960    NOZERO = 0
         GO TO 990

  970    K = TR
         IF (R.LT.N) K = MAX0(K,MINPRM-NZEROK+1)
         IF (ITER.GE.MINPRM .AND. NOZERO.GE.K) GO TO 980
         GO TO 990
C        THEN
C         **TERMINATION BY THE RECURSIVE TEST.
  980    NOCOEF = ITER - NOZERO
         IER = 0
         RETURN

  990    IF (MODE.EQ.1 .AND. ITER.EQ.MINPRM) GO TO 1000
         GO TO 1010
C        THEN
C         **TERMINATION AFTER MINPRM PRIMES.
 1000    NOCOEF = ITER - NOZERO
         IER = 1
         RETURN

 1010 CONTINUE
C   **TERMINATION AFTER MAXPRM PRIMES.
      NOCOEF = MAXPRM - NOZERO
      IER = 1
      RETURN
      END
C     SUBROUTINE VSTORE(A,JA,N,V)
C
C****************************************************************
C                                                               *
C THIS SUBROUTINE COPIES THE VECTOR OF N ELEMENTS               *
C STORED IN THE ARRAY A WITH INCREMENT JA INTO THE VECTOR V     *
C                                                               *
C****************************************************************
C
      SUBROUTINE VSTORE(A,JA,N,V)
      INTEGER A(1),V(1)
      J1 = 1
      DO 10 I = 1,N
         V(I) = A(J1)
         J1 = J1 + JA
   10 CONTINUE
      RETURN
      END
C     INTEGER FUNCTION SCAPRO(A,JA,B,JB,N,P)
C
C****************************************************************
C                                                               *
C THIS FUNCTION COMPUTES THE SCALAR PRODUCT MODULO P            *
C OF TWO INTEGER VECTORS OF N ELEMENTS WHICH ARE STORED         *
C IN THE ARRAYS A AND B WITH INCREMENTS JA AND JB,RESPECTIVELY. *
C                                                               *
C****************************************************************
C
      INTEGER FUNCTION SCAPRO(A,JA,B,JB,N,P)
      INTEGER A(1),B(1),P
      SCAPRO = 0
      J1 = 1
      J2 = 1
      DO 10 I = 1,N
         SCAPRO = MOD(SCAPRO+A(J1)*B(J2),P)
         J1 = J1 + JA
         J2 = J2 + JB
   10 CONTINUE
      RETURN
      END
C     SUBROUTINE SYMSOL(A,NA,N,B,NB,I1,K1,K2,V,CHG,DMOD,P,IER)
C
C****************************************************************
C                                                               *
C THIS SUBROUTINE SOLVES THE SYSTEM OF EQUATIONS AX=B MODULO P  *
C BY COMPUTING                                                  *
C    A**(-1)*B MOD P                                            *
C WHICH IS STORED IN B. ADDITIONALLY                            *
C    DMOD = DETERMINANT(A) MOD P                                *
C IS COMPUTED.                                                  *
C THE MATRIX A OF COEFFICIENTS IS SYMMETRIC (ONLY THE UPPER     *
C TRIANGLE IS NEEDED) AND A SQUARE-ROOT-FREE CHOLESKY           *
C DECOMPOSITION IS USED TO FIND THE SOLUTION. THE PIVOTS ARE    *
C TAKEN FROM THE MAIN DIAGONAL. IF NO NONZERO PIVOT IS          *
C AVAILABLE THE SUBROUTINE RETURNS WITH THE ERROR CODE IER=1.   *
C TO MAKE THE SUBROUTINE UNIVERSALLY APPLICABLE B IS ALLOWED    *
C TO BE ANY CONNECTED PART OF A TWO-DIMENSIONAL ARRAY.          *
C                                                               *
C      A  IS AN INTEGER MATRIX OF DIMENSION NA BY N WHICH       *
C (INPUT) CONTAINS THE MATRIX OF COEFFICIENTS IN THE FIRST      *
C         N ROWS.                                               *
C     NA  IS THE SIZE OF THE FIRST DIMENSION OF A.              *
C (INPUT)                                                       *
C      N  IS THE NUMBER OF EQUATIONS AND THE NUMBER OF          *
C (INPUT) UNKNOWNS IN THE SYSTEM.                               *
C      B  IS AN INTEGER MATRIX OF DIMENSION NB BY K2. IN THE    *
C (INPUT) ROWS I1 TO I1+N-1 AND COLUMNS K1 TO K2 IT CONTAINS    *
C         FIRST THE RIGHT-HAND SIDES B AND FINALLY              *
C         A**(-1)*B MOD P.                                      *
C     NB  IS THE SIZE OF THE FIRST DIMENSION OF B.              *
C (INPUT)                                                       *
C     I1  IS THE NUMBER OF THE FIRST ROW USED IN B.             *
C (INPUT)                                                       *
C     K1  IS THE NUMBER OF THE FIRST COLUMN USED IN B.          *
C (INPUT)                                                       *
C     K2  IS THE NUMBER OF THE LAST COLUMN USED IN B.           *
C (INPUT)                                                       *
C      V  IS A VECTOR OF DIMENSION N USED FOR TEMPORARY         *
C  (TEMP) STORAGE DURING THE BACK SUBSTITUTION.                 *
C    CHG  IS A VECTOR OF DIMENSION N USED FOR SAVING            *
C  (TEMP) SYMMETRIC PERMUTATIONS IN A.                          *
C   DMOD  CONTAINS THE RESIDUUM OF DETERMINANT(A) MODULO P.     *
C(OUTPUT)                                                       *
C      P  IS THE PRIME WHICH IS USED AS MODUL.                  *
C (INPUT)                                                       *
C    IER  IS AN ERROR CODE WHICH IS                             *
C(OUTPUT) 0  IF THE SOLUTION IS COMPUTED.                       *
C         1  IF NOT N NONZERO PIVOTS IN THE MAIN DIAGONAL       *
C            COULD BE FOUND ; NO SOLUTION IS COMPUTED.          *
C         2  IF THE INPUT PARAMETERS ARE INCORRECT.             *
C                                                               *
C                  *****WARNING*****                            *
C                                                               *
C THIS SUBROUTINE ASSUMES THAT THE ELEMENTS OF A AND B ARE      *
C REDUCED MODULO P.                                             *
C                                                               *
C****************************************************************
C
      SUBROUTINE SYMSOL(A,NA,N,B,NB,I1,K1,K2,V,CHG,DMOD,P,IER)
      INTEGER A(NA,N),B(NB,K2),V(N),CHG(N),DMOD,P,PIV
      IER = 2
      I2 = I1 + N - 1
      IF (N.LT.1 .OR. N.GT.NA .OR. I1.LT.1 .OR. I2.GT.NB .OR. K1.LT.
     .    1 .OR. K1.GT.K2) RETURN
      IER = 1
      DMOD = 1
      DO 10 I = 1,N
         CHG(I) = I
   10 CONTINUE
      I = 1
C
C     REPEAT ... I-TH STEP OF THE CHOLESKY DECOMPOSITION OF A
   20 IP1 = I + 1
      IF (A(I,I).EQ.0) GO TO 30
      GO TO 150
C        THEN
C         **LOOK FOR A NONZERO PIVOT IN THE MAIN DIAGONAL.
C
   30 IF (I.EQ.N) RETURN
      DO 40 J = IP1,N
         IF (A(J,J).NE.0) GO TO 50
   40 CONTINUE
C              THEN
C                 TAKE A(J,J) AS PIVOT
C
C           IF NO NONZERO PIVOT WAS FOUND
C           THEN
      RETURN
C           ELSE
C            **INTERCHANGE THE I-TH AND J-TH ROW AND COLUMN OF A.
C
   50 IF (I.GT.1) GO TO 60
      GO TO 80
C              THEN
   60 IM1 = I - 1
      DO 70 K = 1,IM1
         KTEMP = A(K,J)
         A(K,J) = A(K,I)
         A(K,I) = KTEMP
   70 CONTINUE
   80 IF (IP1.LT.J) GO TO 90
      GO TO 110
C              THEN
   90 JM1 = J - 1
      DO 100 K = IP1,JM1
         KTEMP = A(K,J)
         A(K,J) = A(I,K)
         A(I,K) = KTEMP
  100 CONTINUE
  110 IF (J.LT.N) GO TO 120
      GO TO 140
C              THEN
  120 JP1 = J + 1
      DO 130 K = JP1,N
         KTEMP = A(J,K)
         A(J,K) = A(I,K)
         A(I,K) = KTEMP
  130 CONTINUE
  140 KTEMP = A(I,I)
      A(I,I) = A(J,J)
      A(J,J) = KTEMP
      KTEMP = CHG(I)
      CHG(I) = CHG(J)
      CHG(J) = KTEMP
C
  150 PIV = A(I,I)
      DMOD = MOD(DMOD*PIV,P)
      PIV = INV(PIV,P)
      A(I,I) = PIV
      IF (IP1.LE.N) GO TO 160
      GO TO 190
C        THEN
  160 DO 180 J = IP1,N
         KTEMP = MOD(A(I,J)*PIV,P)
         DO 170 K = J,N
            A(J,K) = MOD(A(J,K)-KTEMP*A(I,K),P)
  170    CONTINUE
         A(I,J) = KTEMP
  180 CONTINUE
      I = IP1
      GO TO 20
C        ELSE
C     <-----EXIT
C     CONTINUE
C
C   **SOLVE THE SYSTEM AX=B MOD P USING THE DECOMPOSITION OF A.
C
  190 DO 280 K = K1,K2
         DO 230 I = 1,N
            IB = CHG(I) + I1 - 1
            KTEMP = B(IB,K)
            IF (I.GT.1) GO TO 200
            GO TO 220
C           THEN
  200       IM1 = I - 1
            DO 210 J = 1,IM1
               KTEMP = MOD(KTEMP-A(J,I)*V(J),P)
  210       CONTINUE
  220       V(I) = KTEMP
  230    CONTINUE
         DO 270 I = 1,N
            II = N - I + 1
            KTEMP = MOD(V(II)*A(II,II),P)
            IF (II.LT.N) GO TO 240
            GO TO 260
C           THEN
  240       IP1 = II + 1
            DO 250 J = IP1,N
               KTEMP = MOD(KTEMP-A(II,J)*V(J),P)
  250       CONTINUE
  260       IIB = CHG(II) + I1 - 1
            B(IIB,K) = KTEMP
            V(II) = KTEMP
  270    CONTINUE
  280 CONTINUE
      IER = 0
      RETURN
      END
C     FUNCTION INV(B,P)
C
C****************************************************************
C                                                               *
C EUCLID'S EXTENDED ALGORITHM IS USED TO FIND                   *
C THE INVERSE, INV, MODULO P OF B.                              *
C                                                               *
C                  *****WARNING*****                            *
C                                                               *
C THIS FUNCTION ASSUMES THAT THE GREATEST COMMON DIVISOR        *
C OF B AND P IS 1.                                              *
C                                                               *
C****************************************************************
C
      FUNCTION INV(B,P)
      INTEGER B,P,R0,X0,R1,X1,R2,X2,Q
      R0 = P
      X0 = 0
      R1 = B
      X1 = 1
C     REPEAT
   10 IF (IABS(R1).EQ.1) GO TO 30
      GO TO 20
C        THEN
C     <-----EXIT
   20 Q = R0/R1
      R2 = R0 - Q*R1
      X2 = X0 - Q*X1
      R0 = R1
      X0 = X1
      R1 = R2
      X1 = X2
      GO TO 10
C     CONTINUE
   30 INV = MOD(X1,P)
      IF (R1.LT.0) INV = -INV
      RETURN
      END
C     SUBROUTINE MIXRAD(G,ITER,RES,TEST)
C
C****************************************************************
C                                                               *
C GIVEN THE RESIDUUM RES MODULO P AND THE FIRST ITER-1          *
C COEFFICIENTS IN THE MIXED-RADIX REPRESENTATION                *
C G(1)+G(2)*AP(1)+G(3)*AP(1)*AP(2)+... OF AN INTEGER NUMBER     *
C (WHERE THE AP(I) ARE THE ACTUALLY USED PRIMES)                *
C THIS SUBROUTINE COMPUTES THE ITER-TH COEFFICIENT G(ITER).     *
C                                                               *
C ADDITIONAL PARAMETERS:                                        *
C                                                               *
C     IP  IS AN INTEGER SUCH THAT IP*AP(1)*...*AP(ITER-1) = 1   *
C (INPUT)                                              MODULO P.*
C     P2  IS EQUAL TO (P+1)/2.                                  *
C (INPUT)                                                       *
C   TEST  IS NOT CHANGED IF G(ITER) = 0 WAS COMPUTED;           *
C(OUTPUT) OTHERWISE TEST IS SET .FALSE..                        *
C                                                               *
C****************************************************************
C
      SUBROUTINE MIXRAD(G,ITER,RES,TEST)
      INTEGER G(ITER),RES,P,P2,AP(100)
      LOGICAL TEST
      COMMON /ACTPR/P,IP,P2,AP
      IF (ITER.EQ.1) GO TO 10
      GO TO 20
C     THEN
   10 NTEMP = MOD(RES,P)
      GO TO 60
C     ELSE
   20 ITER1 = ITER - 1
      NTEMP = G(ITER1)
      IF (ITER.NE.2) GO TO 30
      GO TO 50
C        THEN
   30 ITER2 = ITER - 2
      DO 40 I = 1,ITER2
         I1 = ITER1 - I
         NTEMP = MOD(NTEMP*AP(I1)+G(I1),P)
   40 CONTINUE
   50 NTEMP = MOD(IP*MOD(RES-NTEMP,P),P)
   60 IF (NTEMP.NE.0) TEST = .FALSE.
      G(ITER) = NTEMP - (NTEMP/P2)*P
      RETURN
      END
C     SUBROUTINE FRADIX(C,N,A,LMAX,L,B,IER)
C
C****************************************************************
C                                                               *
C GIVEN THE MIXED-RADIX INTEGER                                 *
C                                                               *
C   C(1)+C(2)*AP(1)+ ... +C(N)*AP(1)*...*AP(N-1),               *
C                                                               *
C WHERE AP(I) ARE THE ACTUALLY USED PRIMES AND                  *
C ABS(C(I)) < (P(I)+1)/2, I=1,...,N,                            *
C THIS SUBROUTINE COMPUTES (FOR SOME L <= LMAX) THE COEFFICIENTS*
C A(1),A(2),...,A(L) IN ITS FIXED-RADIX REPRESENTATION:         *
C                                                               *
C   A(1)*B**(L-1)+A(2)*B**(L-2)+...+A(L),                       *
C                                                               *
C WHERE ABS(A(I)) < B AND ALL A(I) HAVE THE SAME SIGN.          *
C                                                               *
C LMAX MUST BE GIVEN SUCH THAT                                  *
C                                                               *
C   LMAX >= (N*LOG(AP(N)) - LOG 2) / LOG B.                     *
C                                                               *
C IER IS AN ERROR CODE WHICH IS 1 IF THE DIMENSION PARAMETERS   *
C N,LMAX ARE INCORRECT, AND 0 OTHERWISE.                        *
C                                                               *
C                     ****** WARNING ******                     *
C                                                               *
C THIS SUBROUTINE ASSUMES THAT FOR ALL I                        *
C                                                               *
C    AP(I)*B                                                    *
C                                                               *
C DOES NOT OVERFLOW A COMPUTER WORD.                            *
C                                                               *
C****************************************************************
C
      SUBROUTINE FRADIX(C,N,A,LMAX,L,B,IER)
      INTEGER C(N),A(LMAX),AP(100),Q,QTEMP,PP,B,P2
      COMMON /ACTPR/PP,IP,P2,AP
      IER = 1
      IF (N.LT.1 .OR. LMAX.LT.1) RETURN
      L = 1
      A(1) = 0
      DO 40 I = 1,N
C        **AT THIS STAGE A(1)+A(2)*B+...+A(L)*B**(L-1) IS THE
C        **FIXED-RADIX REPRESENTATION OF C(N-I+2)+C(N-I+3)*
C        **AP(N-I+2)+...+C(N)*AP(N-I+2)*...*AP(N-1).
         NI1 = N - I + 1
         PP = AP(NI1)
         Q = C(NI1)
C        **COMPUTE THE FIRST L COEFFICIENTS OF THE
C        **FIXED-RADIX REPRESENTATION OF C(N-I+1)+
C        **AP(N-I+1)*(A(1)+A(2)*B+...+A(L)*B**(L-1)).
         DO 10 J = 1,L
            QTEMP = A(J)*PP + Q
            Q = QTEMP/B
            A(J) = QTEMP - Q*B
   10    CONTINUE
C        REPEAT ... CONVERT A(1)+A(2)*B+...+A(L)*B**(L-1)
C               ... +Q*B**L TO FIXED-RADIX FORM.
   20    IF (Q.NE.0) GO TO 30
         GO TO 40
C           THEN
   30    L = L + 1
         IF (L.GT.LMAX) RETURN
         QTEMP = Q/B
         A(L) = Q - QTEMP*B
         Q = QTEMP
         GO TO 20
C        CONTINUE
   40 CONTINUE
      IF (L.GE.2) GO TO 50
      GO TO 70
C     THEN ... REORDER THE COEFFICIENTS.
   50 L2 = L/2
      DO 60 I = 1,L2
         LI1 = L - I + 1
         QTEMP = A(I)
         A(I) = A(LI1)
         A(LI1) = QTEMP
   60 CONTINUE
   70 IER = 0
      RETURN
      END
      BLOCK DATA
C
C****************************************************************
C                                                               *
C                  *****WARNING*****                            *
C                                                               *
C ARRAY NAMES ARE USED IN THE DATA STATEMENTS BELOW TO          *
C SPECIFY VALUES FOR THE ARRAYS KPRIME AND IPRIME. IN SOME      *
C INSTALLATION IT MAY BE NECESSARY TO EXPLICITLY LIST THE       *
C COMPONENTS OF THE ARRAYS, NAMELY, KPRIME(1), KPRIME(2),...,   *
C KPRIME(100) AND IPRIME(1), IPRIME(2),...,IPRIME(100).         *
C                                                               *
C****************************************************************
C
      COMMON /PRIMEB/KPRIME(100),IPRIME(100)
      DATA KPRIME/45233,45247,45259,45263,45281,45289,45293,45307,45317,
     .     45319,45329,45337,45341,45343,45361,45377,45389,45403,45413,
     .     45427,45433,45439,45481,45491,45497,45503,45523,45533,45541,
     .     45553,45557,45569,45587,45589,45599,45613,45631,45641,45659,
     .     45667,45673,45677,45691,45697,45707,45737,45751,45757,45763,
     .     45767,45779,45817,45821,45823,45827,45833,45841,45853,45863,
     .     45869,45887,45893,45943,45949,45953,45959,45971,45979,45989,
     .     46021,46027,46049,46051,46061,46073,46091,46093,46099,46103,
     .     46133,46141,46147,46153,46171,46181,46183,46187,46199,46219,
     .     46229,46237,46261,46271,46273,46279,46301,46307,46309,46327,
     .     46337/
      DATA IPRIME/00000,42015,28577,01108,29342,16641,10405,19447,26685,
     .     39525,14116,12753,32178,01043,08857,27911,15049,07079,33425,
     .     00804,23175,23886,44779,41942,10171,16606,10638,17371,27195,
     .     35827,42639,01829,24658,09023,37958,30638,06339,41270,40538,
     .     10157,11783,00457,32947,42170,17910,33474,20017,25086,36508,
     .     37444,35543,06993,10326,16328,26765,42083,37223,30711,09408,
     .     06635,38421,11397,32683,17333,34245,15748,35735,23492,19302,
     .     20076,45620,44978,09864,14832,16092,19457,24045,44950,32872,
     .     24309,15726,43057,37766,14046,41826,19946,41363,23967,39791,
     .     29237,18085,12952,36850,02213,30023,34871,42667,40410,32615,
     .     46136/
      END