File: SIMQ.FT of Tape: Various/ETH/eth11-1
(Source file text)
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 C