C
C     ..................................................................
C
C        SUBROUTINE SIMQ
C
C        PURPOSE
C           OBTAIN SOLUTION OF A SET OF SIMULTANEOUS LINEAR EQUATIONS,
C           AX=B
C
C        USAGE
C           CALL SIMQ(A,B,N,KS)
C
C        DESCRIPTION OF PARAMETERS
C           A - MATRIX OF COEFFICIENTS STORED COLUMNWISE.  THESE ARE
C               DESTROYED IN THE COMPUTATION.  THE SIZE OF MATRIX A IS
C               N BY N.
C           B - VECTOR OF ORIGINAL CONSTANTS (LENGTH N). THESE ARE
C               REPLACED BY FINAL SOLUTION VALUES, VECTOR X.
C           N - NUMBER OF EQUATIONS AND VARIABLES. N MUST BE .GT. ONE.
C           KS - OUTPUT DIGIT
C                0 FOR A NORMAL SOLUTION
C                1 FOR A SINGULAR SET OF EQUATIONS
C
C        REMARKS
C           MATRIX A MUST BE GENERAL.
C           IF MATRIX IS SINGULAR , SOLUTION VALUES ARE MEANINGLESS.
C           AN ALTERNATIVE SOLUTION MAY BE OBTAINED BY USING MATRIX
C           INVERSION (MINV) AND MATRIX PRODUCT (GMPRD).
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           METHOD OF SOLUTION IS BY ELIMINATION USING LARGEST PIVOTAL
C           DIVISOR. EACH STAGE OF ELIMINATION CONSISTS OF INTERCHANGING
C           ROWS WHEN NECESSARY TO AVOID DIVISION BY ZERO OR SMALL
C           ELEMENTS.
C           THE FORWARD SOLUTION TO OBTAIN VARIABLE N IS DONE IN
C           N STAGES. THE BACK SOLUTION FOR THE OTHER VARIABLES IS
C           CALCULATED BY SUCCESSIVE SUBSTITUTIONS. FINAL SOLUTION
C           VALUES ARE DEVELOPED IN VECTOR B, WITH VARIABLE 1 IN B(1),
C           VARIABLE 2 IN B(2),........, VARIABLE N IN B(N).
C           IF NO PIVOT CAN BE FOUND EXCEEDING A TOLERANCE OF 0.0,
C           THE MATRIX IS CONSIDERED SINGULAR AND KS IS SET TO 1. THIS
C           TOLERANCE CAN BE MODIFIED BY REPLACING THE FIRST STATEMENT.
C
C     ..................................................................
C
      SUBROUTINE SIMQ(A,B,N,KS)
      DIMENSION A(1),B(1)
C
C        FORWARD SOLUTION
C
      TOL=0.0
      KS=0
      JJ=-N
      DO 65 J=1,N
      JY=J+1
      JJ=JJ+N+1
      BIGA=0
      IT=JJ-J
      DO 30 I=J,N
C
C        SEARCH FOR MAXIMUM COEFFICIENT IN COLUMN
C
      IJ=IT+I
      IF(ABS(BIGA)-ABS(A(IJ))) 20,30,30
   20 BIGA=A(IJ)
      IMAX=I
   30 CONTINUE
C
C        TEST FOR PIVOT LESS THAN TOLERANCE (SINGULAR MATRIX)
C
      IF(ABS(BIGA)-TOL) 35,35,40
   35 KS=1
      RETURN
C
C        INTERCHANGE ROWS IF NECESSARY
C
   40 I1=J+N*(J-2)
      IT=IMAX-J
      DO 50 K=J,N
      I1=I1+N
      I2=I1+IT
      SAVE=A(I1)
      A(I1)=A(I2)
      A(I2)=SAVE
C
C        DIVIDE EQUATION BY LEADING COEFFICIENT
C
   50 A(I1)=A(I1)/BIGA
      SAVE=B(IMAX)
      B(IMAX)=B(J)
      B(J)=SAVE/BIGA
C
C        ELIMINATE NEXT VARIABLE
C
      IF(J-N) 55,70,55
   55 IQS=N*(J-1)
      DO 65 IX=JY,N
      IXJ=IQS+IX
      IT=J-IX
      DO 60 JX=JY,N
      IXJX=N*(JX-1)+IX
      JJX=IXJX+IT
   60 A(IXJX)=A(IXJX)-(A(IXJ)*A(JJX))
   65 B(IX)=B(IX)-(B(J)*A(IXJ))
C
C        BACK SOLUTION
C
   70 NY=N-1
      IT=N*N
      DO 80 J=1,NY
      IA=IT-J
      IB=N-J
      IC=N
      DO 80 K=1,J
      B(IB)=B(IB)-A(IA)*B(IC)
      IA=IA-N
   80 IC=IC-1
      RETURN
      END
