C ACM ALGORITHM 575 C C PERMUTATIONS FOR A ZERO-FREE DIAGONAL C C BY I.S. DUFF C C ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, SEPTEMBER, 1981 C C TEST DECK FOR MC21A ... SI/ 1 C THESE RUNS ON RANDOM MATRICES CAUSE ALL POSSIBLE STATEMENTS IN 2 C MC21A TO BE EXECUTED. 3 INTEGER IP(21) 4 C INTEGER*2 ICN(1000),IPERM(20),LENR(20),IW1(80) I/ 5 INTEGER ICN(1000),IPERM(20),LENR(20),IW1(80) 6 LOGICAL SING 7 COMMON /RANMA/ SING 8 LICN=1000 9 C RUN ON RANDOM MATRICES OF ORDERS 1 THROUGH 20. 10 DO 100 N=1,20 11 NP1=N+1 12 NOV4=MAX0(N/4,1) 13 10 CALL FA01BS(NOV4,L) 14 KNUM=L*N 15 C KNUM IS REQUESTED NUMBER OF NON-ZEROS IN RANDOM MATRIX. 16 IF (KNUM.GT.LICN) GO TO 10 17 C IF SING IS .FALSE., MATRIX IS GUARANTEED STRUCTURALLY NON-SINGULAR. 18 SING=(N/2)*2.EQ.N 19 C CALL TO SUBROUTINE TO GENERATE RANDOM MATRIX. 20 CALL RANMAT(N,N,ICN,LICN,IP,NP1,KNUM,IW1) 21 C KNUM IS NOW ACTUAL NUMBER OF NON-ZEROS IN RANDOM MATRIX. 22 IF (KNUM.GT.LICN) GO TO 10 23 WRITE(6,20) N,KNUM,SING 24 20 FORMAT(//5H N = ,I2,6H NZ = ,I4,8H SING = ,L1) 25 C SET UP ARRAY OF ROW LENGTHS. 26 DO 30 I=1,N 27 30 LENR(I)=IP(I+1)-IP(I) 28 C CALL TO MC21A. 29 CALL MC21A(N,ICN,LICN,IP,LENR,IPERM,NUMNZ,IW1) 30 C 31 C TESTING TO SEE IF THERE ARE NUMNZ NON-ZEROS ON THE DIAGONAL 32 C OF THE PERMUTED MATRIX. 33 NUM=0 34 DO 50 I=1,N 35 IOLD=IPERM(I) 36 J1=IP(IOLD) 37 J2=J1+LENR(IOLD)-1 38 IF (J2.LT.J1) GO TO 50 39 DO 40 JJ=J1,J2 40 J=ICN(JJ) 41 IF (J.NE.I) GO TO 40 42 NUM=NUM+1 43 GO TO 50 44 40 CONTINUE 45 50 CONTINUE 46 IF (NUM.NE.NUMNZ) WRITE(6,60) NUMNZ,NUM 47 60 FORMAT(28H FAILURE IN MC21A, NUMNZ IS ,I2,11H INSTEAD OF,I3) 48 C 49 100 CONTINUE 50 C STOP END C SUBROUTINE TO GENERATE RANDOM MATRIX. 54 SUBROUTINE RANMAT(M,N,ICN,LICN,IPTR,NNNP1,KNUM,IW) 55 INTEGER IPTR(NNNP1) C INTEGER*2 ICN(LICN),IW(N) I/ INTEGER ICN(LICN),IW(N) LOGICAL SING COMMON /RANMA/ SING INUM=(KNUM/N)*2 IF (INUM.GT.N-1) INUM=N-1 MATNUM=1 C EACH PASS THROUGH THIS LOOP GENERATES A ROW OF THE MATRIX. DO 50 J=1,M IPTR(J)=MATNUM IF (SING) GO TO 10 IF (J.GT.N) GO TO 10 ICN(MATNUM)=J MATNUM=MATNUM+1 10 IF (N.EQ.1) GO TO 50 DO 20 I=1,N 20 IW(I)=0 IF (.NOT.SING) IW(J)=1 CALL FA01BS(INUM,LROW) LROW=LROW-1 IF (LROW.EQ.0) GO TO 50 C LROW OFF-DIAGONAL NON-ZEROS IN ROW J OF THE MATRIX. DO 40 II=1,LROW 30 CALL FA01BS(N,I) IF (IW(I).EQ.1) GO TO 30 IW(I)=1 ICN(MATNUM)=I MATNUM=MATNUM+1 40 CONTINUE 50 CONTINUE MP1=M+1 DO 60 I=MP1,NNNP1 60 IPTR(I)=MATNUM KNUM=MATNUM-1 RETURN END C RANDOM NUMBER GENERATOR. 93 FUNCTION FA01AS(I) 94 COMMON/FA01ES/G DOUBLE PRECISION G G=DMOD(G* 9228907.D0,4294967296.D0) IF(I.GE.0)FA01AS=G/4294967296.D0 IF(I.LT.0)FA01AS=2.D0*G/4294967296.D0-1.D0 RETURN END SUBROUTINE FA01BS(MAX,NRAND) 102 NRAND=INT(FA01AS(1)*FLOAT(MAX))+1 RETURN END BLOCK DATA 106 COMMON/FA01ES/G 107 DOUBLE PRECISION G 108 DATA G/1431655765.D0/ 109 END 110 C S IS THE STANDARD FORTRAN VERSION SI/ 1 C I IS THE IBM VERSION 2 C 3 SUBROUTINE MC21A(N,ICN,LICN,IP,LENR,IPERM,NUMNZ,IW) 4 C C DESCRIPTION OF PARAMETERS. C INPUT VARIABLES N,ICN,LICN,IP,LENR. C OUTPUT VARIABLES IPERM,NUMNZ. C C N ORDER OF MATRIX. C ICN ARRAY CONTAINING THE COLUMN INDICES OF THE NON-ZEROS. THOSE C BELONGING TO A SINGLE ROW MUST BE CONTIGUOUS BUT THE ORDERING C OF COLUMN INDICES WITHIN EACH ROW IS UNIMPORTANT AND WASTED C SPACE BETWEEN ROWS IS PERMITTED. C LICN LENGTH OF ARRAY ICN. C IP IP(I), I=1,2,...N, IS THE POSITION IN ARRAY ICN OF THE FIRST C COLUMN INDEX OF A NON-ZERO IN ROW I. C LENR LENR(I) IS THE NUMBER OF NON-ZEROS IN ROW I, I=1,2,..N. C IPERM CONTAINS PERMUTATION TO MAKE DIAGONAL HAVE THE SMALLEST C NUMBER OF ZEROS ON IT. ELEMENTS (IPERM(I),I) I=1, ... N ARE C NON-ZERO AT THE END OF THE ALGORITHM UNLESS MATRIX C IS STRUCTURALLY SINGULAR. IN THIS CASE, (IPERM(I),I) WILL C BE ZERO FOR N-NUMNZ ENTRIES. C NUMNZ NUMBER OF NON-ZEROS ON DIAGONAL OF PERMUTED MATRIX. C IW WORK ARRAY .. SEE LATER COMMENTS. C INTEGER IP(N) C INTEGER*2 ICN(LICN),LENR(N),IPERM(N),IW(N,4) I/ INTEGER ICN(LICN),LENR(N),IPERM(N),IW(N,4) CALL MC21B(N,ICN,LICN,IP,LENR,IPERM,NUMNZ,IW(1,1),IW(1,2),IW(1,3), 1IW(1,4)) RETURN END C/ 34 SUBROUTINE MC21B(N,ICN,LICN,IP,LENR,IPERM,NUMNZ,PR,ARP,CV,OUT) 35 INTEGER IP(N) C C DIVISION OF WORK ARRAY IS NOW DESCRIBED. C C PR(I) IS THE PREVIOUS ROW TO I IN THE DEPTH FIRST SEARCH. C ARP(I) IS ONE LESS THAN THE NUMBER OF NON-ZEROS IN ROW I C WHICH HAVE NOT BEEN SCANNED WHEN LOOKING FOR A CHEAP ASSIGNMENT. C CV(I) IS THE MOST RECENT ROW EXTENSION AT WHICH COLUMN I C WAS VISITED. C OUT(I) IS ONE LESS THAN THE NUMBER OF NON-ZEROS IN ROW I C WHICH HAVE NOT BEEN SCANNED DURING ONE PASS THROUGH THE C MAIN LOOP. C C INTEGER*2 ICN(LICN),LENR(N),IPERM(N),PR(N),CV(N), I/ C 1ARP(N),OUT(N) I/ INTEGER ICN(LICN),LENR(N),IPERM(N),PR(N),CV(N), 1ARP(N),OUT(N) C C INITIALIZATION OF ARRAYS. DO 10 I=1,N ARP(I)=LENR(I)-1 CV(I)=0 10 IPERM(I)=0 NUMNZ=0 C C C MAIN LOOP. C EACH PASS ROUND THIS LOOP EITHER RESULTS IN A NEW ASSIGNMENT C OR GIVES A ROW WITH NO ASSIGNMENT. DO 130 JORD=1,N J=JORD PR(J)=-1 DO 100 K=1,JORD C LOOK FOR A CHEAP ASSIGNMENT IN1=ARP(J) IF (IN1.LT.0) GO TO 60 IN2=IP(J)+LENR(J)-1 IN1=IN2-IN1 DO 50 II=IN1,IN2 I=ICN(II) IF (IPERM(I).EQ.0) GO TO 110 50 CONTINUE C NO CHEAP ASSIGNMENT IN ROW. ARP(J)=-1 C BEGIN LOOKING FOR ASSIGNMENT CHAIN STARTING WITH ROW J. 60 OUT(J)=LENR(J)-1 C INNER LOOP. EXTENDS CHAIN BY ONE OR BACKTRACKS. DO 90 KK=1,JORD IN1=OUT(J) IF (IN1.LT.0) GO TO 80 IN2=IP(J)+LENR(J)-1 IN1=IN2-IN1 C FORWARD SCAN. DO 70 II=IN1,IN2 I=ICN(II) IF (CV(I).EQ.JORD) GO TO 70 C COLUMN I HAS NOT YET BEEN ACCESSED DURING THIS PASS. J1=J J=IPERM(I) CV(I)=JORD PR(J)=J1 OUT(J1)=IN2-II-1 GO TO 100 70 CONTINUE C C BACKTRACKING STEP. 80 J=PR(J) IF (J.EQ.-1) GO TO 130 90 CONTINUE C 100 CONTINUE C C NEW ASSIGNMENT IS MADE. 110 IPERM(I)=J ARP(J)=IN2-II-1 NUMNZ=NUMNZ+1 DO 120 K=1,JORD J=PR(J) IF (J.EQ.-1) GO TO 130 II=IP(J)+LENR(J)-OUT(J)-2 I=ICN(II) IPERM(I)=J 120 CONTINUE C 130 CONTINUE C C IF MATRIX IS STRUCTURALLY SINGULAR, WE NOW COMPLETE THE C PERMUTATION IPERM. IF (NUMNZ.EQ.N) GO TO 500 DO 140 I=1,N 140 ARP(I)=0 K=0 DO 160 I=1,N IF (IPERM(I).NE.0) GO TO 150 K=K+1 OUT(K)=I GO TO 160 150 J=IPERM(I) ARP(J)=I 160 CONTINUE K=0 DO 170 I=1,N IF (ARP(I).NE.0) GO TO 170 K=K+1 IOUTK=OUT(K) IPERM(IOUTK)=I 170 CONTINUE 500 RETURN END