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*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