C ACM ALGORITHM 576 C C A FORTRAN PROGRAM FOR SOLVING AX = B C C BY I. BARRODALE AND G.F. STUART C C ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, SEPTEMBER, 1981 C C THIS MAIN PROGRAM CONTAINS TEN TEST PROGRAMS. MAN 10 C THE FIRST FOUR ILLUSTRATE THE USE OF THE ALGORITHM, MAN 20 C THE LAST SIX ARE SIMPLE TESTS TO EXECUTE ALL MAN 30 C STATEMENTS OF THE SUBROUTINES. MAN 40 COMPLEX Z(35) MAN 50 REAL R(32) MAN 60 REAL AS(32,32), BS(32), XP(32), S(32), DX(32) MAN 70 REAL A(32,32), B(32), X(32), RZ(32) MAN 80 REAL CHAR(20) MAN 90 INTEGER P(32), Q(32) MAN 100 KIN = 5 MAN 110 KOUT = 6 MAN 120 C TEST 1 FROM BARRODALE AND STEWART J. INST. MATHS APPLICS MAN 130 C (1977) -V19-, P39-47. MAN 140 TOL = 1.E-6 MAN 150 EPS = 1.0 MAN 160 TOL = 0 MAN 170 N = 3 MAN 180 ND = 32 MAN 190 READ (KIN,99999) (CHAR(I),I=1,20) MAN 200 WRITE (KOUT,99998) (CHAR(I),I=1,20) MAN 210 DO 10 I=1,N MAN 220 READ (KIN,99997) (A(I,J),J=1,N), B(I) MAN 230 WRITE (KOUT,99996) (A(I,J),J=1,N), B(I) MAN 240 10 CONTINUE MAN 250 DO 30 I=1,N MAN 260 BS(I) = B(I) MAN 270 DO 20 J=1,N MAN 280 AS(I,J) = A(I,J) MAN 290 20 CONTINUE MAN 300 30 CONTINUE MAN 310 CALL MODGE(N, ND, A, B, EPS, TOL, P, Q, RZ, X, K, MAN 320 * IOC, IOP) MAN 330 CALL OUTPT1(K, IOC, IOP, X, N) MAN 340 CALL RESID1(AS, BS, X, N, ND) MAN 350 ERR = 1.E-7 MAN 360 ITER = 50 MAN 370 CALL REFINE(N, ND, AS, BS, A, P, Q, K, ERR, R, S, XP, MAN 380 * DX, X, ITER, DIGITS) MAN 390 CALL OUTPT2(ITER, DIGITS, X, N) MAN 400 CALL RESID2(AS, BS, X, N, ND) MAN 410 C TEST 2 A POLYNOMIAL INTERPOLATION OF AN MAN 420 C EXPONENTIAL FORM. MAN 430 TOL = 1.E-6 MAN 440 EPS = 0.0 MAN 450 N = 5 MAN 460 ND = 32 MAN 470 DO 50 I=1,N MAN 480 XX = I*0.25 MAN 490 DO 40 KK=1,N MAN 500 A(I,KK) = XX**KK MAN 510 40 CONTINUE MAN 520 B(I) = XX*EXP(-XX) MAN 530 50 CONTINUE MAN 540 DO 70 I=1,N MAN 550 BS(I) = B(I) MAN 560 DO 60 J=1,N MAN 570 AS(I,J) = A(I,J) MAN 580 60 CONTINUE MAN 590 70 CONTINUE MAN 600 CALL MODGE(N, ND, A, B, EPS, TOL, P, Q, RZ, X, K, MAN 610 * IOC, IOP) MAN 620 CALL OUTPT1(K, IOC, IOP, X, N) MAN 630 CALL RESID1(AS, BS, X, N, ND) MAN 640 ERR = 1.E-7 MAN 650 ITER = 50 MAN 660 CALL REFINE(N, ND, AS, BS, A, P, Q, K, ERR, R, S, XP, MAN 670 * DX, X, ITER, DIGITS) MAN 680 CALL OUTPT2(ITER, DIGITS, X, N) MAN 690 CALL RESID2(AS, BS, X, N, ND) MAN 700 C TEST 3 A RANK DEFICIENT INTERPOLATION OF AN MAN 710 C EXPONENTIAL FORM. MAN 720 TOL = 1.E-6 MAN 730 EPS = 1.D-6 MAN 740 N = 9 MAN 750 ND = 32 MAN 760 DO 80 I=1,N MAN 770 XX = (I-1)/8.0 MAN 780 A(I,1) = 1.0 MAN 790 A(I,2) = XX MAN 800 A(I,3) = (XX-1.0) MAN 810 A(I,4) = XX**2 MAN 820 A(I,5) = (XX**2-XX) MAN 830 A(I,6) = XX**3 MAN 840 A(I,7) = (XX**3-XX**2) MAN 850 A(I,8) = XX**4 MAN 860 A(I,9) = (XX**4-XX**3) MAN 870 B(I) = EXP(XX) MAN 880 80 CONTINUE MAN 890 DO 100 I=1,N MAN 900 BS(I) = B(I) MAN 910 DO 90 J=1,N MAN 920 AS(I,J) = A(I,J) MAN 930 90 CONTINUE MAN 940 100 CONTINUE MAN 950 CALL MODGE(N, ND, A, B, EPS, TOL, P, Q, RZ, X, K, * IOC, IOP) CALL OUTPT1(K, IOC, IOP, X, N) CALL RESID1(AS, BS, X, N, ND) ERR = 1.E-7 ITER = 50 CALL REFINE(N, ND, AS, BS, A, P, Q, K, ERR, R, S, XP, * DX, X, ITER, DIGITS) CALL OUTPT2(ITER, DIGITS, X, N) CALL RESID2(AS, BS, X, N, ND) C TEST 4 A COMPLEX INTERPOLATION PROBLEM ARISING FROM A C DIRICHLET PROBLEM. TOL = 1.E-6 EPS = 0.0 EPS = 0.05 N = 32 ND = 32 READ (KIN,99999) (CHAR(I),I=1,20) WRITE (KOUT,99998) (CHAR(I),I=1,20) READ (KIN,99995) (Z(I),I=1,N) WRITE (KOUT,99994) (Z(I),I=1,N) DO 120 I=1,N A(I,1) = 1.0 DO 110 J=1,15 A(I,J*2) = REAL(Z(I)**J) A(I,J*2+1) = AIMAG(Z(I)**J) 110 CONTINUE A(I,32) = REAL(Z(I)**16) SS = REAL(Z(I)) T = AIMAG(Z(I)) B(I) = EXP(SS**2-SS*T+2.0*T**2) 120 CONTINUE DO 140 I=1,N BS(I) = B(I) DO 130 J=1,N AS(I,J) = A(I,J) 130 CONTINUE 140 CONTINUE CALL MODGE(N, ND, A, B, EPS, TOL, P, Q, RZ, X, K, * IOC, IOP) CALL OUTPT1(K, IOC, IOP, X, N) CALL RESID1(AS, BS, X, N, ND) ERR = 1.E-7 ITER = 50 CALL REFINE(N, ND, AS, BS, A, P, Q, K, ERR, R, S, XP, * DX, X, ITER, DIGITS) CALL OUTPT2(ITER, DIGITS, X, N) CALL RESID2(AS, BS, X, N, ND) C C THE FOLLOWING ARE THE SIX SIMPLE TESTS. C C SIMPLE TEST 1 (THE A MATRIX IS ZERO). TOL = 1.E-6 EPS = 0.0 ND = 32 N = 1 A(1,1) = 0.D0 B(1) = 1.D0 AS(1,1) = A(1,1) CALL MODGE(N, ND, A, B, EPS, TOL, P, Q, RZ, X, K, * IOC, IOP) CALL OUTPT1(K, IOC, IOP, X, N) C SIMPLE TEST 2 (A PROBLEM WHICH NEEDS NO SOLUTION OR C REFINEMENT). TOL = 1.E-6 EPS = 0.0 ND = 32 N = 1 A(1,1) = 1.D0 B(1) = 1.D0 AS(1,1) = A(1,1) CALL MODGE(N, ND, A, B, EPS, TOL, P, Q, RZ, X, K, * IOC, IOP) CALL OUTPT1(K, IOC, IOP, X, N) CALL RESID1(AS, BS, X, N, ND) ERR = 1.E-7 ITER = 1 CALL REFINE(N, ND, AS, BS, A, P, Q, K, ERR, R, S, XP, * DX, X, ITER, DIGITS) CALL OUTPT2(ITER, DIGITS, X, N) CALL RESID2(AS, BS, X, N, ND) C SIMPLE TEST 3 (A PROBLEM OF SIZE ZERO). TOL = 1.E-6 EPS = 0.0 ND = 32 N = 0 A(1,1) = 0.D0 B(1) = 1.D0 AS(1,1) = A(1,1) CALL MODGE(N, ND, A, B, EPS, TOL, P, Q, RZ, X, K, * IOC, IOP) CALL OUTPT1(K, IOC, IOP, X, N) C SIMPLE TEST 4 (A TEST OF SMALL NUMBERS). TOL = 1.E-6 EPS = 1.E-8 N = 2 DO 160 I=1,N DO 150 J=1,N A(I,J) = 0.5*TOL 150 CONTINUE B(I) = EPS 160 CONTINUE CALL MODGE(N, ND, A, B, EPS, TOL, P, Q, RZ, X, K, * IOC, IOP) CALL OUTPT1(K, IOC, IOP, X, N) C SIMPLE TEST 5 -PART 1- (TEST SMALL PIVOTS). EPS = 1.E-6 TOLER = 1.E-5 A(1,1) = 1.E0 A(2,1) = 0.E00 A(1,2) = 0.E00 A(2,2) = 1.E-7 B(1) = 1.E0 B(2) = 1.E0 + 1.E-7 N = 2 CALL MODGE(N, ND, A, B, EPS, TOLER, P, Q, RZ, X, K, * IOC, IOP) CALL OUTPT1(K, IOC, IOP, X, N) CALL REFINE(N, ND, AS, BS, A, P, Q, K, ERR, R, S, XP, * DX, X, ITER, DIGITS) C SIMPLE TEST 5 -PART 2- (REFINE A SOLUTION OF SIZE ZERO). K = 0 CALL REFINE(N, ND, AS, BS, A, P, Q, K, ERR, R, S, XP, * DX, X, ITER, DIGITS) C SIMPLE TEST 6 (TEST 3 ABOVE USED TO TEST SMALL PIVOT). TOL = 1.E-7 EPS = 1.E-6 N = 9 ND = 32 DO 170 I=1,N XX = (I-1)/8.0 A(I,1) = 1.0 A(I,2) = XX A(I,3) = (XX-1.0) A(I,4) = XX**2 A(I,5) = (XX**2-XX) A(I,6) = XX**3 A(I,7) = (XX**3-XX**2) A(I,8) = XX**4 A(I,9) = (XX**4-XX**3) B(I) = EXP(XX) 170 CONTINUE DO 190 I=1,N BS(I) = B(I) DO 180 J=1,N AS(I,J) = A(I,J) 180 CONTINUE 190 CONTINUE CALL MODGE(N, ND, A, B, EPS, TOL, P, Q, RZ, X, K, * IOC, IOP) STOP 99999 FORMAT (20A1) 99998 FORMAT (12H1 INPUT FOR , 20A1) 99997 FORMAT (4F5.0) 99996 FORMAT (2H , 4E20.7) 99995 FORMAT (2F5.3) 99994 FORMAT (3H Z=, 2E20.7) END SUBROUTINE MODGE(N, NDIM, A, B, EPS, TOLER, P, Q, Z, MOD 10 * X, K, OCODE, OPIVOT) C THIS SUBROUTINE SOLVES ANY N-BY-N SYSTEM OF LINEAR C EQUATIONS AX=B, USING A MODIFICATION OF GAUSSIAN C ELIMINATION. THE MODIFICATION OCCURS IN THE CHOICE OF C PIVOTS: AT EACH STAGE THE PIVOT ROW IS DETERMINED BY THE C ABSOLUTELY LARGEST RESIDUAL (RIGHT HAND SIDE), AND THEN C THE PIVOT IS CHOSEN AS THE ABSOLUTELY LARGEST COEFFICIENT C IN THE PIVOT ROW. THE PIVOT IS POSITIONED ON THE PRINCIPAL C DIAGONAL BY (IF NECESSARY) A ROW INTERCHANGE AND A COLUMN C INTERCHANGE. C C IF AFTER THE K-TH STAGE THE LARGEST RESIDUAL IS TOO SMALL C (.LT.EPS), THE ELIMINATION IS TERMINATED. BACK C SUBSTITUTION IS THEN USED TO OBTAIN AN APPROXIMATE C SOLUTION WHICH SATISFIES THE FIRST K EQUATIONS. C C IF THE (K+1)-TH PIVOT IS TOO SMALL (.LE.TOLER), THE C ELIMINATION IS TERMINATED. BACK SUBSTITUTION IS THEN USED C TO OBTAIN AN APPROXIMATE SOLUTION WHICH SATISFIES THE C FIRST K EQUATIONS. C C DESCRIPTION OF PARAMETERS: C N AN INTEGER DENOTING THE NUMBER OF EQUATIONS C (N.GE.1). C NDIM AN INTEGER PARAMETER FOR ADJUSTABLE DIMENSIONS C (NDIM.GE.N). C A A REAL ARRAY OF NDIM ROWS AND AT LEAST N C COLUMNS WHICH ON ENTRY MUST CONTAIN THE COEFFICIENT C MATRIX A IN ITS FIRST N ROWS AND N COLUMNS. C ON EXIT THE FIRST K ROWS AND K COLUMNS (K.LE.N) OF C A CONTAIN THE LU DECOMPOSITION OF THE K EQUATIONS C SATISFIED. C B A REAL ARRAY OF DIMENSION AT LEAST N WHICH C ON ENTRY MUST CONTAIN IN POSITIONS 1 TO N THE RIGHT C HAND SIDE VECTOR B. C ON EXIT THE POSITIONS K+1 TO N CONTAIN THE RESIDUALS C OF THE UNSATISFIED EQUATIONS. C EPS A REAL INPUT TOLERANCE. IF THE (ABSOLUTELY) LARGEST C RESIDUAL AT ANY STAGE OF THE ELIMINATION IS LESS C THAN EPS THEN THE ELIMINATION IS TERMINATED AND C BACK SUBSTITUTION COMMENCES. PREMATURE TERMINATION C OF THIS TYPE CAN BE AVOIDED BY SETTING EPS TO ZERO. C TOLER A REAL POSITIVE INPUT TOLERANCE. IF THE ABSOLUTE C VALUE OF THE PIVOT AT ANY STAGE OF THE ELIMINATION C IS LESS THAN OR EQUAL TO TOLER, THEN COMPLETE C PIVOTING IS USED FOR THIS STAGE. IF THE ABSOLUTE C VALUE OF THE NEW PIVOT IS LESS THAN OR EQUAL TO C TOLER, THE REMAINING EQUATIONS ARE CONSIDERED C TO BE LINEARLY DEPENDENT, AND SO THE ELIMINATION C IS TERMINATED AND BACK SUBSTITUTION COMMENCES. C TOLER SHOULD NORMALLY BE SET TO APPROXIMATELY C 10**(-D+1), WHERE D REPRESENTS THE NUMBER OF C DECIMAL DIGITS OF ACCURACY AVAILABLE. HOWEVER, C PROBLEMS WHICH ARE NOT 'REASONABLY SCALED' C BEFOREHAND COULD REQUIRE A MUCH LARGER OR C SMALLER VALUE FOR TOLER THAN IS RECOMMENDED C HERE. C P AN INTEGER ARRAY OF DIMENSION AT LEAST N USED TO C RECORD COLUMN INTERCHANGES OF THE ARRAY A. C ON EXIT P(I) IS THE INDEX OF THE ORIGINAL POSITION C OF COLUMN I. C Q AN INTEGER ARRAY OF DIMENSION AT LEAST N USED TO C RECORD ROW INTERCHANGES OF THE ARRAYS A AND B. C ON EXIT Q(I) IS THE INDEX OF THE ORIGINAL POSITION C OF ROW I. C Z A REAL WORK SPACE ARRAY OF DIMENSION AT LEAST N. C X A REAL OUTPUT ARRAY OF DIMENSION AT LEAST N C WHICH CONTAINS THE SOLUTION (IN CORRECT ORDER) OF C THE N LINEAR EQUATIONS IN POSITIONS 1 TO N. C IF THE SOLUTION IS OBTAINED USING LESS THAN N C EQUATIONS, THEN THE REMAINING COMPONENTS OF X ARE C SET TO ZERO. C K AN INTEGER WHICH ON OUTPUT DENOTES THE NUMBER OF C EQUATIONS SATISFIED (K.LE.N). C OCODE AN INTEGER EXIT CODE WITH VALUES: C 0- EITHER N.LE.0 OR (N.EQ.1 AND A(1,1).EQ.0). C 1- THE SYSTEM AX=B IS SOLVED. C 2- THE ELIMINATION IS TERMINATED BY USE OF C EPS (SMALL RESIDUAL). C 3- THE ELIMINATION IS TERMINATED BY USE OF C TOLER (SMALL PIVOT). C OPIVOT AN INTEGER EXIT CODE WITH VALUES: C 1- COMPLETE PIVOTING IS NOT EMPLOYED. C 2- COMPLETE PIVOTING IS EMPLOYED, HENCE SMALL C RESIDUALS (.LT.EPS) MAY HAVE BEEN ALLOWED. C C AN IMPLEMENTATION NOTE: C IF YOUR FORTRAN COMPILER PERMITS A SINGLE COLUMN OF A TWO C DIMENSIONAL ARRAY TO BE PASSED TO A ONE DIMENSIONAL ARRAY C THROUGH A SUBROUTINE CALL, CONSIDERABLE SAVINGS IN C EXECUTION TIME MAY BE ACHIEVED BY USING SUBROUTINE C SAXPY (SEE A.C.M. TRANSACTIONS ON MATHEMATICAL C SOFTWARE VOL.5 NUM.3 1979 PP.308-323) C (SEE STATEMENT NUMBERS 160, 180, AND 220 BELOW.) DIMENSION A(NDIM,1), B(1), X(1), Z(1) INTEGER P(1), Q(1), OCODE, OPIVOT DATA IONE /1/ IF (N.GT.0) GO TO 10 OCODE = 0 GO TO 290 10 OPIVOT = 1 IF (N.GT.1) GO TO 30 IF (A(1,1).EQ.0.0) GO TO 20 X(1) = B(1)/A(1,1) P(1) = 1 Q(1) = 1 K = 1 OCODE = 1 GO TO 290 20 OCODE = 0 GO TO 290 30 NM1 = N - 1 C INITIALIZE COLUMN PERMUTATION AND ROW PERMUTION VECTORS. DO 40 I=1,N P(I) = I Q(I) = I 40 CONTINUE C THE FOLLOWING DO LOOP IS THE MAJOR LOOP IN THE ELIMINATION C PROCEDURE. DO 190 K=1,NM1 KP1 = K + 1 C DETERMINE THE NEXT EQUATION TO BE SATISFIED, I.E. THE C EQUATION WITH THE LARGEST RESIDUAL. MROW = K RESMAX = ABS(B(MROW)) DO 50 I=KP1,N IF (RESMAX.GE.ABS(B(I))) GO TO 50 MROW = I RESMAX = ABS(B(MROW)) 50 CONTINUE C IF THE LARGEST RESIDUAL IS LESS THAN EPS THEN START THE C BACK SUBSTITUTION. IF (RESMAX.GE.EPS) GO TO 60 M = K OCODE = 2 GO TO 210 C DETERMINE THE PIVOT COLUMN, I.E. THE COLUMN WHICH CONTAINS C THE ABSOLUTELY LARGEST COEFFICIENT IN THE PIVOT ROW C (THE EQUATION TO BE SATISFIED). 60 MCOL = K AMAX = ABS(A(MROW,MCOL)) DO 70 J=KP1,N IF (AMAX.GE.ABS(A(MROW,J))) GO TO 70 MCOL = J AMAX = ABS(A(MROW,J)) 70 CONTINUE C IF THE ABSOLUTE VALUE OF THE PIVOT IS GREATER THAN TOLER C THEN CONTINUE WITH THE ELIMINATION, OTHERWISE PERFORM C COMPLETE PIVOTING ON THE ELEMENTS OF THE UNSATISFIED C EQUATIONS TO DETERMINE A NEW PIVOT. IF (AMAX.GT.TOLER) GO TO 100 AMAX = 0. DO 90 J=K,N DO 80 I=K,N IF (AMAX.GT.ABS(A(I,J))) GO TO 80 MROW = I MCOL = J AMAX = ABS(A(MROW,MCOL)) 80 CONTINUE 90 CONTINUE OPIVOT = 2 C IF THE ABSOLUTE VALUE OF THE PIVOT RESULTING FROM COMPLETE C PIVOTING IS LESS THAN OR EQUAL TO TOLER THEN START THE C BACK SUBSTITUTION. IF (AMAX.GT.TOLER) GO TO 100 M = K OCODE = 3 GO TO 210 C IF THE PIVOT ROW IS NOT THE K-TH ROW THEN PERFORM C THE NECESSARY PERMUTATIONS ON THE ELEMENTS OF A AND B C TO MAKE IT SO AND RECORD THESE PERMUTATIONS IN Q. 100 IF (MROW.EQ.K) GO TO 120 DO 110 J=1,N T = A(K,J) A(K,J) = A(MROW,J) A(MROW,J) = T 110 CONTINUE T = B(K) B(K) = B(MROW) B(MROW) = T M = Q(K) Q(K) = Q(MROW) Q(MROW) = M C IF THE PIVOT COLUMN IS NOT THE K-TH COLUMN THEN PERFORM C THE NECESSARY PERMUTATIONS ON A TO MAKE IT SO AND RECORD C THESE PERMUTATIONS IN P. 120 IF (MCOL.EQ.K) GO TO 140 DO 130 I=1,N T = A(I,K) A(I,K) = A(I,MCOL) A(I,MCOL) = T 130 CONTINUE M = P(K) P(K) = P(MCOL) P(MCOL) = M C RECORD THE MULTIPLIERS FOR THE K-TH ELIMINATION. 140 T = -A(K,K) DO 150 I=KP1,N A(I,K) = A(I,K)/T 150 CONTINUE C PERFORM THE K-TH ELIMINATION STEP IN A COLUMN-ORIENTED C MANNER. DO 170 J=KP1,N T = A(K,J) C IF POSSIBLE, USE THE SUBROUTINE SAXPY AND REPLACE THE C FOLLOWING THREE STATEMENTS WITH: C INDX=N-KP1+1 C 160 CALL SAXPY (INDX,T,A(KP1,K),IONE,A(KP1,J),IONE) DO 160 I=KP1,N A(I,J) = A(I,J) + A(I,K)*T 160 CONTINUE 170 CONTINUE T = B(K) C IF POSSIBLE, USE THE SUBROUTINE SAXPY AND REPLACE THE C FOLLOWING THREE STATEMENTS WITH: C INDX=N-KP1+1 C 180 CALL SAXPY (INDX,T,A(KP1,K),IONE,B(KP1),IONE) DO 180 I=KP1,N B(I) = B(I) + A(I,K)*T 180 CONTINUE 190 CONTINUE M = N OCODE = 2 IF (ABS(B(N)).LT.EPS) GO TO 200 M = N + 1 OCODE = 1 200 IF (ABS(A(N,N)).GT.TOLER) GO TO 210 M = N OCODE = 3 210 K = M - 1 IF (K.EQ.0) GO TO 250 C PERFORM BACK SUBSTITUTION ON THE K EQUATIONS C SATISFIED, IN A COLUMN-ORIENTED MANNER. IF (K.EQ.1) GO TO 240 KM1 = K - 1 DO 230 JJ=1,KM1 JM1 = K - JJ J = JM1 + 1 Z(J) = B(J)/A(J,J) T = -Z(J) C IF POSSIBLE, USE THE SUBROUTINE SAXPY AND REPLACE THE C FOLLOWING THREE STATEMENTS WITH: C INDX=JM1 C 220 CALL SAXPY (INDX,T,A(1,J),IONE,B,IONE) DO 220 I=1,JM1 B(I) = B(I) + A(I,J)*T 220 CONTINUE 230 CONTINUE 240 Z(1) = B(1)/A(1,1) IF (K.EQ.N) GO TO 270 250 KP1 = K + 1 C ZERO THE VARIABLES WHICH WERE NOT USED. DO 260 I=KP1,N Z(I) = 0. 260 CONTINUE C REORDER THE VARIABLES TO CORRESPOND TO THE ORDER OF THE C INPUT EQUATIONS. 270 DO 280 I=1,N M = P(I) X(M) = Z(I) 280 CONTINUE 290 RETURN END SUBROUTINE SAXPY(N, SA, SX, INCX, SY, INCY) SAX 10 C C OVERWRITE SINGLE PRECISION SY WITH SINGLE PRECISION SA*SX +SY. C FOR I = 0 TO N-1, REPLACE SY(LY+I*INCY) WITH SA*SX(LX+I*INCX) + C SY(LY+I*INCY), WHERE LX = 1 IF INCX .GE. 0, ELSE LX = (-INCX)*N, C AND LY IS DEFINED IN A SIMILAR WAY USING INCY. C REAL SX(N), SY(N), SA IF (N.LE.0 .OR. SA.EQ.0.E0) RETURN IF (INCX.EQ.INCY) IF (INCX-1) 10, 30, 70 10 CONTINUE C C CODE FOR NONEQUAL OR NONPOSITIVE INCREMENTS. C IX = 1 IY = 1 IF (INCX.LT.0) IX = (-N+1)*INCX + 1 IF (INCY.LT.0) IY = (-N+1)*INCY + 1 DO 20 I=1,N SY(IY) = SY(IY) + SA*SX(IX) IX = IX + INCX IY = IY + INCY 20 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 4. C 30 M = MOD(N,4) IF (M.EQ.0) GO TO 50 DO 40 I=1,M SY(I) = SY(I) + SA*SX(I) 40 CONTINUE IF (N.LT.4) RETURN 50 MP1 = M + 1 DO 60 I=MP1,N,4 SY(I) = SY(I) + SA*SX(I) SY(I+1) = SY(I+1) + SA*SX(I+1) SY(I+2) = SY(I+2) + SA*SX(I+2) SY(I+3) = SY(I+3) + SA*SX(I+3) 60 CONTINUE RETURN C C CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS. C 70 CONTINUE NS = N*INCX DO 80 I=1,NS,INCX SY(I) = SA*SX(I) + SY(I) 80 CONTINUE RETURN END SUBROUTINE RESID1(A, B, X, N, ND) RES 10 DOUBLE PRECISION SUM, DBLE DIMENSION A(ND,1), B(1), X(1) KOUT = 6 WRITE (KOUT,99999) DO 20 I=1,N SUM = 0.D0 DO 10 J=1,N SUM = SUM + DBLE(A(I,J))*DBLE(X(J)) 10 CONTINUE RESIDU = B(I) - SUM WRITE (KOUT,99998) RESIDU 20 CONTINUE RETURN 99999 FORMAT (15H RESIDUALS ARE ) 99998 FORMAT (E18.7) END SUBROUTINE RESID2(A, B, X, N, ND) RES 10 DOUBLE PRECISION SUM, DBLE DIMENSION A(ND,1), B(1), X(1) KOUT = 6 WRITE (KOUT,99999) DO 20 I=1,N SUM = 0.D0 DO 10 J=1,N SUM = SUM + DBLE(A(I,J))*DBLE(X(J)) 10 CONTINUE RESIDU = B(I) - SUM WRITE (KOUT,99998) RESIDU 20 CONTINUE RETURN 99999 FORMAT (28H RESIDUALS AFTER REFINE ARE ) 99998 FORMAT (E18.7) END SUBROUTINE REFINE(N, NDIM, A, B, ALU, P, Q, KMAX, REF 10 * ERR, R, S, XNEW, DX, X, ITER, DIGITS) C THIS SUBROUTINE PERFORMS ITERATIVE IMPROVEMENT ON THE C COMPUTED SOLUTION X OBTAINED FROM SUBROUTINE MODGE. C THE ITERATIONS CONTINUE UNTIL EITHER THE CORRECTIONS C TO X ARE NEGLIGIBLE, I.E. CONVERGENCE OCCURS (SEE C DESCRIPTION OF ERR), OR UNTIL COMPLETION OF THE C MAXIMUM NUMBER OF ITERATIONS ALLOWED (ITER). C C THE FOLLOWING TABLE GIVES THE CORRESPONDENCE BETWEEN THE C OUTPUT PARAMETERS OF MODGE AND THE INPUT PARAMETERS OF C REFINE: C NAME ON OUTPUT FROM MODGE: NAME ON INPUT TO REFINE: C N N C NDIM NDIM C A ALU C P P C Q Q C K KMAX C X X C C DESCRIPTION OF PARAMETERS: C N AN INTEGER DENOTING THE NUMBER OF EQUATIONS C (N.GE.1). C NDIM AN INTEGER PARAMETER FOR ADJUSTABLE DIMENSIONS C (NDIM.GE.N). C A A REAL ARRAY OF NDIM ROWS AND AT LEAST N COLUMNS C WHICH ON ENTRY MUST CONTAIN THE COEFFICIENT C MATRIX A IN ITS FIRST N ROWS AND N COLUMNS. C A IS NOT ALTERED BY REFINE. C B A REAL ARRAY OF DIMENSION AT LEAST N WHICH ON C ENTRY MUST CONTAIN IN POSITIONS 1 TO N THE RIGHT C HAND SIDE VECTOR B. B IS NOT ALTERED BY REFINE. C ALU THE REAL OUTPUT ARRAY A OF SUBROUTINE MODGE. ALU IS C NOT ALTERED BY REFINE. C P AN INTEGER ARRAY OF DIMENSION AT LEAST N, USED C TO RECORD COLUMN INTERCHANGES OF THE ARRAY A. C P IS THE OUTPUT ARRAY P OF SUBROUTINE MODGE. C Q AN INTEGER ARRAY OF DIMENSION AT LEAST N, USED C TO RECORD ROW INTERCHANGES OF THE ARRAYS A AND B. C Q IS THE OUTPUT ARRAY Q OF SUBROUTINE MODGE. C KMAX THE INTEGER OUTPUT PARAMETER K OF MODGE. C ERR A REAL INPUT PARAMETER USED TO TEST FOR C CONVERGENCE. ERR SHOULD NOT BE LESS THAN 10**(-D), C WHERE D REPRESENTS THE NUMBER OF DECIMAL DIGITS C OF ACCURACY AVAILABLE. C THE CONVERGENCE CRITERION IS: C EXIT IF THE FOLLOWING RELATIONSHIP IS TRUE FOR C ALL I=1 TO KMAX: C ABS(DX(I)).LE.ABS(ERR*X(I)).OR.X(I).EQ.0.0 . C R A REAL WORK SPACE ARRAY OF DIMENSION AT LEAST N. C (USED BY REFINE TO CALCULATE THE RESIDUALS AT EACH C ITERATION.) C S A REAL WORK SPACE ARRAY OF DIMENSION AT LEAST N. C (USED BY REFINE TO CALCULATE THE CORRECTIONS, DX, C AT EACH ITERATION.) C XNEW A REAL WORK SPACE ARRAY OF THE SAME DIMENSION AS X. C DX A REAL WORK SPACE ARRAY OF DIMENSION AT LEAST N. C (USED BY REFINE TO CALCULATE A CORRECTION TO X C AT EACH ITERATION.) C X THE REAL OUTPUT ARRAY X OF SUBROUTINE MODGE. C ON EXIT X CONTAINS (USUALLY) A BETTER SOLUTION C TO THE LINEAR SYSTEM AX=B. C ITER AN INTEGER DENOTING, ON INPUT, THE MAXIMUM NUMBER C OF ITERATIONS OF ITERATIVE IMPROVEMENT ALLOWED. C ITER SHOULD NOT EXCEED 2*D, WHERE D REPRESENTS C THE NUMBER OF DECIMAL DIGITS OF ACCURACY AVAILABLE. C ON OUTPUT ITER IS THE ACTUAL NUMBER OF ITERATIONS C COMPLETED. C DIGITS A REAL PARAMETER WHICH ON OUTPUT INDICATES THE C APPROXIMATE NUMBER OF SIGNIFICANT DIGITS IN THE C INITIAL SOLUTION (I.E. THE SOLUTION FROM C MODGE). C C AN IMPLEMENTATION NOTE: C IF YOUR FORTRAN COMPILER PERMITS A SINGLE COLUMN OF A TWO C DIMENSIONAL ARRAY TO BE PASSED TO A ONE DIMENSIONAL ARRAY C THROUGH A SUBROUTINE CALL, CONSIDERABLE SAVINGS IN C EXECUTION TIME MAY BE ACHIEVED BY USING SUBROUTINE SAXPY C (SEE ABOVE) WHICH OPERATES ON COLUMN VECTORS. C (SEE ALSO STATEMENT NUMBERS 50 AND 60 BELOW.) DOUBLE PRECISION SUM, DBLE C SUM IS USED IN THE DOUBLE PRECISION COMPUTATION OF THE C RESIDUALS. DIMENSION A(NDIM,1), ALU(NDIM,1), B(1), X(1), * XNEW(1), DX(1), R(1), S(1) INTEGER P(1), Q(1) LOGICAL FIN C FIN IS A LOGICAL VARIABLE WHICH IS USED TO TEST FOR C CONVERGENCE OF ITERATIVE IMPROVEMENT. DATA IDEFLT, EDEFLT /14,1.E-7/ C IDEFLT AND EDEFLT ARE MACHINE DEPENDENT DEFAULT VALUES FOR C ITER AND ERR. THEY ARE APPROPRIATE FOR AN IBM SYSTEM 360/ C 370 AS SHOWN. IN GENERAL THEY SHOULD BE SET, RESPECTIVELY, C TO 2*D AND 10**(-D) WHERE D REPRESENTS THE NUMBER OF C DECIMAL DIGITS OF ACCURACY AVAILABLE. DATA IONE /1/ ITMAX = MIN0(ITER,IDEFLT) EPS = AMAX1(ERR,EDEFLT) DIGITS = -ALOG10(EDEFLT) KMAXP1 = KMAX + 1 KMAXM1 = KMAX - 1 ITER = 0 IF (KMAX.EQ.0) GO TO 190 C FILL THE PERMUTED ARRAY R=BP. DO 10 J=1,N II = Q(J) R(J) = B(II) XNEW(J) = X(J) DX(J) = XNEW(J) 10 CONTINUE C THE FOLLOWING DO LOOP IS THE MAJOR LOOP IN REFINE. DO 160 IT=1,ITMAX C FORM THE RESIDUALS IN DOUBLE PRECISION. DO 30 I=1,N SUM = 0.0D0 DO 20 J=1,N IF (Q(J).EQ.I) JJ = J SUM = SUM + DBLE(A(I,J))*DBLE(DX(J)) 20 CONTINUE R(JJ) = DBLE(R(JJ)) - SUM 30 CONTINUE DO 40 I=1,KMAX S(I) = R(I) 40 CONTINUE IF (KMAX.EQ.1) GO TO 90 C APPLY THE ELIMINATION PROCEDURE TO THE RIGHT HAND SIDES C (I.E. THE RESIDUALS) IN A COLUMN ORIENTED MANNER. DO 60 K=1,KMAXM1 KP1 = K + 1 T = S(K) C IF POSSIBLE, USE THE SUBROUTINE SAXPY AND REPLACE THE C FOLLOWING THREE STATEMENTS WITH: C INDX=KMAX-KP1+1 C 50 CALL SAXPY (INDX,T,ALU(KP1,K),IONE,S(KP1),IONE) DO 50 I=KP1,KMAX S(I) = S(I) + ALU(I,K)*T 50 CONTINUE 60 CONTINUE C FORM THE CORRECTIONS TO THE SOLUTION BY BACK SUBSTITUTION, C IN A COLUMN ORIENTED MANNER. DO 80 JJ=1,KMAXM1 JM1 = KMAX - JJ J = JM1 + 1 DX(J) = S(J)/ALU(J,J) T = -DX(J) C IF POSSIBLE, USE THE SUBROUTINE SAXPY AND REPLACE THE C FOLLOWING THREE STATEMENTS WITH: C INDX=JM1 C 70 CALL SAXPY (INDX,T,ALU(1,J),IONE,S,IONE) DO 70 I=1,JM1 S(I) = S(I) + ALU(I,J)*T 70 CONTINUE 80 CONTINUE 90 DX(1) = S(1)/ALU(1,1) C TEST FOR COMPLETION AND FORM THE CORRECTED SOLUTION. FIN = .TRUE. IF (KMAX.GE.N) GO TO 110 DO 100 I=KMAXP1,N DX(I) = 0.0 100 CONTINUE 110 DO 120 I=1,N II = P(I) S(II) = DX(I) 120 CONTINUE DO 130 I=1,N DX(I) = S(I) T = X(I) IF (T.EQ.0.0) GO TO 130 IF (ABS(DX(I)).GT.ABS(EPS*T)) FIN = .FALSE. XNEW(I) = XNEW(I) + DX(I) 130 CONTINUE C IF (.NOT.FIN) GO TO 140 ITER = IT GO TO 170 140 IF (IT.NE.1) GO TO 160 C ESTIMATE THE NUMBER OF SIGNIFICANT DIGITS IN THE SOLUTION. DO 150 I=1,N IF (X(I).EQ.0.0 .OR. DX(I).EQ.0.0) GO TO 150 DIGITS = AMIN1(DIGITS,ALOG10(ABS(X(I)/DX(I)))) 150 CONTINUE 160 CONTINUE ITER = ITMAX C STORE THE REFINED SOLUTION IN X. 170 DO 180 I=1,N X(I) = XNEW(I) 180 CONTINUE 190 RETURN END SUBROUTINE OUTPT2(ITER, DIGITS, X, N) OUT 10 DIMENSION X(1) KOUT = 6 WRITE (KOUT,99999) WRITE (KOUT,99998) (X(I),I=1,N) WRITE (KOUT,99997) ITER, DIGITS RETURN 99999 FORMAT (40H THE SOLUTION VECTOR X AFTER REFINE IS ) 99998 FORMAT (E18.7) 99997 FORMAT (12H ITERATIONS, I5, 12H DIGITS , E18.8) END SUBROUTINE OUTPT1(K, IOC, IOP, X, N) OUT 10 DIMENSION X(1) KOUT = 6 WRITE (KOUT,99998) IOC, IOP IF (IOC.EQ.0) GO TO 10 WRITE (KOUT,99997) K WRITE (KOUT,99996) WRITE (KOUT,99995) (X(I),I=1,N) GO TO 20 10 WRITE (6,99999) 20 RETURN 99999 FORMAT (13H NOTE OCODE=0, //, 20H NO SOLUTION RETURNE, * 1HD) 99998 FORMAT (9H1 OUTCODE, I5, 15H PIVOTCODE , I5) 99997 FORMAT (38H THE NUMBER OF EQUATIONS SATISIFIED IS, I5) 99996 FORMAT (27H THE SOLUUTION VECTOR X IS ) 99995 FORMAT (E18.7) END TEST 1 DATA READ IN DATA 1.0 1.0 1.0 89.0 DATA 1.0 4.0 16.0209.0 DATA 1.0 5.0 25.0201.0 DATA TEST 4 DATA READ IN DATA 0.0 0.0 DATA 0.0 0.125 DATA 0.0 0.25 DATA 0.0 0.375 DATA 0.0 0.5 DATA 0.0 0.625 DATA 0.0 0.75 DATA 0.0 0.875 DATA 0.0 1.0 DATA 0.1251.0 DATA 0.25 1.0 DATA 0.3751.0 DATA 0.5 1.0 DATA 0.6251.0 DATA 0.75 1.0 DATA 0.8751.0 DATA 1.0 1.0 DATA 1.0 0.875 DATA 1.0 0.75 DATA 1.0 0.625 DATA 1.0 0.5 DATA 1.0 0.375 DATA 1.0 0.25 DATA 1.0 0.125 DATA 1.0 0.0 DATA 0.8750.0 DATA 0.75 0.0 DATA 0.6250.0 DATA 0.5 0.0 DATA 0.3750.0 DATA 0.25 0.0 DATA 0.1250.0 DATA 0.0 0.0 DATA 0.125 DATA 0.25 DATA 0.375 DATA 0.5 DATA 0.625 DATA 0.75 DATA 0.875 DATA 1.0 DATA 1.0 0.125 DATA 1.0 0.25 DATA 1.0 0.375 DATA 1.0 0.50 DATA 1.0 0.625 DATA 1.0 0.75 DATA 1.0 0.875 DATA 1.0 1.0 DATA 0.8751.0 DATA 0.75 1.0 DATA 0.6251.0 DATA 0.5 1.0 DATA 0.3751.0 DATA 0.25 1.0 DATA 0.1251.0 DATA 0.0 1.0 DATA 0.0 0.875 DATA 0.0 0.75 DATA 0.0 0.625 DATA 0.0 0.5 DATA 0.0 0.375 DATA 0.0 0.25 DATA 0.0 0.125 DATA