Time - Space Tradeoff in the Implementation of the Hellerman-Rarick Algorithm by Yonatan Levy and Alexander Meeraus Technical Note No. 10 Research Project 671-58 March 1979 Development Research Center World Bank 1818 H Street, N.W. Washington, D. C. 20433 Preliminary and Confidential: Not for quotations or attributions withov prior clearance from the authors. Views expressed are those of the atchors and do not necessarily reflect those c- the World Bank. 1. Ceneral Descripticn The Hellerman-Rarick (HR) algorithm is a heuristic routine that assigns pivot elements and spikes so as to minimize the number of columns with nonzeros abcve the diagonal. A detailed description of the algorithm is given in Hellerman and Rarick [1]. In this paper we focus our attention on the implementation. We present two alternatives: one that reduces storage space as much as possible, and the other which concentrates on fast execution whle using several auxiliary arrays. The Algorithm The version of the hR algorithm we are dealing with operates on a square irreducible matrix. This matrix is usually a submatrix (bump) of a large sparse matrix. At each step the algorithm stacks a spike or assigns a pivot. Whenever a pivot is assigned, some spikes may be unstacked and assigned to their final positions. The two measures which govern the selection process are the nonzeros count for the rows and the tally function for the columns. The algorith performs as many steps as the order of the matrix. Figure 1 shows the general flow of the algorit-. Every step starts with the colum. selection which is based on th2 maxim-- tally function. The exact rule for breaking ties depends on the implementaticn. Once a column is chosen, the minimum row count is checked. If it is two or more, the selected column is stacked as a spike, and after some updating the step is done. If, on the other hand, the minimum row count is one, the sel2cted column becomes a pivot column. By the method of selection, this col=n is guaranteed to have at least one nonzero in a "count one" row, and this row becomes the pivot row. 3 Entry Initialization Select one column with maximum tally Stacl- the coluMn Assign the column to as a spike. ou count - 1? a pivot position. Update arrays. Update arrays. Any free Exit columns left? you count = 0? Unstacl: a spike Figure 1: General Flo-w-Chart of the Hellerman-Rarick Algorithm If thL column has additional nonzeros in "count one" rows, these rows are used for unstacking spikes. Some updating is necessary before the next step can be carried out. jm2lementation The algorithm was designed for sparse matrices, and such matrices are usually stored in a sparse form, either ro'--wise or column-wise. Hence, only columns (or rows) can -) scanned. For example, if a matrix is stored column-wise, and one would like to update the tally values after a row has been assigned, one could not do it unless a row-wise representation is avail- able. So evidently, additional arrays may reduce the execution time of the algorithm considerably. The various implementations of the algorithm vary in the form the matrix Is stored, in the way the information on tally function values and nonzeros counts is kept, in the rules for breaking ties between maximum tally columns, and, as a result of all these, in the updating that has to be done at the end of each step. The two versions that are offered in this work present two extremes. One, the subroutine .RS, uses the global colu=n-wise represen- tation and the only additional array is needed for storing the row counts. Due to the lack of row-wise arrays the lexical maximum tally selection rule (described in section 4.1 of [2]) is used. This rule uses a weighted tally function which has to be recomputed every step. This selection procedure is relatively slow, but there is a little updating to do. The other version, subroutine HRF, uses a local row-wise representation in addition to the global column-wise arrays. It also keeps updated arrays for row counts, column counts and tally functions. Rows and columns are kept sorted within the permutation arrays, and this requires two additional arrays for pointers. As a result of all that subroutine HRF requires some set-up time, and most of the execution time is devoted to updating. On the other hand, the column selec- tion is very efficient, and over all ERF is much fasttr than HRS, especially for large bumps. ?ermutation Arrays The actual output of the algorithm is a new order of the columns and rows of the bump. Permutation arrays are used to input the initial order which is a result of the block lower triangularization ((21, section 3). The same permutation arrays are updated throughout the algorithm, and in subrou- tine HRF they are used to keep sorted lists as well. On exit, the permutation arrays contain the final order of columns and rows. The arrays are global and reflect the order of all the columns and rows of the matrix. The subroutines are actually operating on a specific set of en:ries in the array. A general discussion on the subject of permutation vectors and matrices is given in section 5 of [2]. Any other array refers to columns or rows by their original indices. Therefore, when columns are exchanged only the column permutation arrays have to be updated. 2. The Subroutines In this section we describe briefly the role of each of the subrou- tines. There exists a definite order in the way these subroutines call ea:h other, as illustrated by the tree in Figure 2. The subroutines are described in the same order a breadth first search rould scan the tree. 6 n Program P42 BUMP ROWREP HRiF FR LSORT L_TALLY REPRC NA:TLY Figure 2: Procedence graph of the subroutines. The root, subroutine ?4F, is a version of subroutine P4X (described in [2], section 4). It was modified to use a work array IW with variable dimension NW. As a result of this m.3dification subroutine BLNP is called instead of SPIKE. We recall that subroutine P4F has two cask. It clusters successive I x 1 blocks to one diagonal block, and it processes the bumps. For each bump it first calls subroutine PARTBL to partition the columns of the bump, then subroutine BUMP is called. Subroutine BUMP is a snort linking subroutine that assigns the execution of the Hellerman-Rarick algorithm to one of two routines. In section I we have described the ti=e-space tradeoff and the two subroutines 4RF and HRS. The work array IW and its size NW are transferred to BUMP as arguments. Subroutine HRF needs seven arrays on cop of those required by lRS. The role of subroutine BUMP is to check whether these arrays could fit into the space ot the array 1W. If they fit, subroutine ROWREP is called first to generate the row-wise representation of the bump. Then BUMP calculates the pointers to the work array 1*, and calls subroutine hIRF. It is assumed that the column permutarton is loaded into the forward positions of IW before BUMP is called. This i crucial if subroutine ERF is used. Whenever the array 1W is not big enough, subroutine ?RS is called. In chis case, the work array is not used at all. Subroutine HRS is a modified -.ersion of subroutine SPIY ([21, section 4.1). The modification made i: slightly more efficient, anc similar in struc,ure to subroutine HRF. The zeneral structure is described in section 1. Since elements can be scanned by ccl=n only, the selection rule for che column is the weighted tally functica. - 8 - Subroutine MLAXTLY chooses the column vith the a2zi=un tally as follows: For each free column j let mk(j) = number of nonzeros which lie in rows with curreat count exactly k, Mk(j) = number of nonzeros in rows with current count of at least k, and define 8 T(j) - 1025 (j) + 1019-2k (j) Q) k=1 where MIN is the current minimum row count. The free column which has the maximal T(j) is selected. This selection procedure is time consuming, especially when the bu=p is large. On the other hand, the only updating that has to be done bet.een steps is to reduce the row counts, and to find the minimum row coc=c after a pivot has been assigned. Subrouti.. HRF Laplements closely the algortch2 as it is described in [1]. It uses both column--ise and row-wise represencation, all four permuta- tion arrays (HRS uses only three), aLd several other at=iTiary arrays in order to execute faster. The permutation arrays are used for storing :he final assignments of the columns and rows, for stacking the spikes, and also for keeping the free rows and columns sorted. The 'ows are initially sarted by subroutine LSORT according to their nonzeros count. Whenever a colu-2 is stacked Cr assigned, the row counts are reduced. This is doae by subroutine REDRC -hizh also keeps the free rows ordered. The columns are scerced ` subroutine TALLY according to their tally function values. "!ienever it is called TALLY recomputes the tally values, and simultaneously creates the sorted list - 9- of free columns. (The tally function is NMIN(j), i.e. the number of non- zeros a column has in minimum count rows). Once the tally values are computed, one column is selected in the following manner: (i) If a unique column has the maximum tally value, this column is selected. (ii) If all the free rows have the same row count, the last column in the list is chosen arbitrarily. (Since the list is sorted, the last column is known to have maximum tally.) (iii) If neither (i) nor (ii) holds, and the maximum tally value is one, subroutine MAXTLY is used to break the tie among the maximum tally columns (and not among all free columns as in subroutine HRS). (iv) If the maaximum tally value is more than one, the column with with maximum nonzeros count is selected. After one column is selected, it either becomes a spike or a pivot column. From this point the program follows the flow chart in Figure 1. The updating is done by subroutines REDRC and TALLY. Subroutine REDRC is called exactly once every step, so it is called as many times as the order of the bump. Subroutine TALLY is not called when there is a unique pivot row. In this case the pivot cclumn is uniquely determined, and the selection procedure is not used. However, except for this special case subroutine TALLY is called once every step. Subroutine REDRC scans only one column, and subroutine TALLY scans only the rows with minimum row count. Subroutine MAXTLY is called only when several columns are tied with maximum tally equals one, and then cnly these columns are scanned. So on the average, subroutine HRF is performing ..ess computations every step than subroutine RRS. The latter calls subroutine MAXTLY every step, and it scans all the free columns. Hence, tne difference between the two subroutines is especially big in the first steps and for big bumps. Exact timings are given in the following section. 3. Performance The whole sequence of subroutines has been tested on several random matrices. One test has used Subroutine HRS for the HR rou'ine, while the other has used Subroutine HRF. The results are summarized in Table 1. From this table one can learn what is the additional storage required by Subroutine HRF, and what is the exact reduction in xecution time. The execution tie is the elapsed time between calling and returning from Subourtine P4F. It is clear from the table that The difference between the two implementations is considerable in matrices with large bumps. For small bumps the efficiency of HRF is balanced by the set-up and updating time. Table 1: Storage Requirements and Execution Times Total No. Columns Total No. Size,Nonzeros,Spikes Storage Execution times - P4F* Matrix Size of Nonzeros Bumps Lu Bumps of Spikes** of Largest Bump** for IIRF Using IIRF Using HRS 1 500 5411 2 375 137 342,3326,133 6013 .117 1.146 2 500 5356 5 276 45 114,557,23 2788 .047 .148 3 1000 10725 2 754 263 682,6467,252 14622 .335 4.483 4 1000 6245 3 436 53 300,1182,36 5785 .106 .547 5 1000 6327 8 213 21 65,179,3 4312 .034 .064 6 1000 10870 8 551 90 226,1154,42 5609 .126 .597 7 2000 12172 1 1429 290 1429,8410,290 192'1 1.529 13.142 8 2000 12522 7 1090 147 438,2210,79 11089 .334 2.120 9 2000 21877 9 624 57 212,725,20 9154 .143 .555 10 2000 21659 16 266 29 128,443,10 8702 .058 .113 F- * La sOcos8, on Lt CDC Cyber 176. ** Spikes counts are based on subroutine IIRF. - 12 - References [1 Hellerman, E., and Rarick, D., Reinversions with the preassigned Aivot routine. Math. ?rog. 1, pp. 195-216 (1971). (21 Meeraus A., and Levy, Y., Programs for Structure Analysis of Sparse Matrices. Tech. Rep. No. 8, Res. Proj. 671-58, DRC, World Bank (1978). Appendix 3: Listings of the Programs A.1 SUBROUTINE P4F (M,RI,CS,CNTO,RNTO,ROTN,BLKS,IB, P4F 1 HSPIKE,CE,RC,IW,NW) P4F C P4F C P4F IS A VERSION OF THE P4 ROUTINE OF HELLERMAN AND RARICK. P4F C IT ACCEPTS AS INPUT A MATRIX IN BLOCK LOWER TRIANGULAR FORM P4F C ( AS CREATED BY SUBROUTINE BLKTR ). P4F GROUPS SUCCESSIVE P4F C ONE-ELEMENT-BLOCKS INTO ONE TRIANGULAR BLOCK, AND OPERATES P4F C ON THE BUMPS ( BY THE SUBROUTINE BUMP ) TO REDUCE THE NUMBER P4F C CF SPIKES. THE OUTPUT INCLUDES THE FINAL PERMUTATIONS, AND P4F C ALL THE INFORMATION ON THE BLOCKS STRUCTURE AND THE SPIKES. P4F C P4F INTEGER RI(1),CS(1),CNTO(1),RNTO(1),ROTN(1),BLKS(1) P4F INTEGER HSPIKE(1),CE(1),RC(1),IW(NW) P4F C P4F C VARIABLES ( THE SUBSCRIPTS ARE THE REQUIRED DIMENSIONS) P4F C P4F C INPUT ARGUMENTS P4F C M ORDER OF THE MATRIX P4F C RI( ) COLUMNWISE ROW INDICES LIST. REQUIRED LENGTH IS CS(M+1). P4F C CS(M+1) ?OINTERS TO BEGINNING OF COLUMNS IN RI( ) P4F C P4F C INPUT/OUT?UT ARGUMENTS P4F C CNTO(M) INVERSE COLUMN PERMUTATION. CNTO(I) IS THE ORIGINAL P4F C INDEX OF THE COLUMN IN THE I'TH POSITION. P4F C RNTO(M) INVERSE ROW PERMUTATION P4F C ROTN(M) THE ROW PERMUTATION P4F C BLKS(M) POINTERS TO STARTING POSITIONS OF BLOCKS. P4F C AT ENTRY BLKS(M-IB+J) IS THE START POSITION OF THE J'TH P4F C IRREDUCIBLE BLOCK. UPON RETURN BLKS(J) IS THE START OF P4F C THE J'TH BLOCK ( EITHER TRIANGULAR OR IRREDUCIBLE ). P4F C IB THE NUMBER OF BLOCKS. THIS NUMBER IS LIKELY TO DECREASE P4F C DURING THE PROGRAM. P4F C P4F C OUTPUT ARGUMENTS P4F C HSPIKE(M) IF HSPIKE(J)=O, THEN COLUMN J IS NOT A SPIKE. P4F C OTHERWISE HSPIKE(J) IS THE HEIGHT OF THE SPIKE. OR - P4F C IF J IS THE FIRST COLUMN OF THE BUMP- HSPIKE(J) IS THE P4F C NUM3ER OF SPIKES IN THE BUMP. P4F C CE(M) ?OINTERS TO THE END OF BUMP IN RI( ) (FOR BUMP COLUMNS) P4F C RC(M) WORK ARRAY USED FOR ROWS NONZEROS COUNT P)F C IW(NW) WORK ARRAY. USED IN SUBROUTINE BUMP. IF IW( ) IS P4F C BIG ENOUGH, THE BUMP IS PROCESSED BY A FASTER ROUTINE, P4F C P4F U SUBROUTNES REQUIRED P).F C PARTBL TO PARTITION A SPECIFIC RANGE OF COLUMNS P4F C BUMP TO FIND THE SPIKES STRUCTURE WITHIN A BUMP Pl-F C P4F C REMARKS P4F C 1. ALL PERMUTATION ARRAYS ARE KEPT UPDATED THROUGHOUT THE ?ROGRAM. P4F C 2. ALL ARRAYS, EXCEPT THE INVERSE PERMUIATION VECTORS, ARE PlF C INDEXED SY THE ORIGINAL INDICES. P1LF C P4F C P4F C SET INITIAL VALUES P4F C P4F LENTR=O P4F MPl=M+1 P4F A.2 IS=MP1-IB P4F BLKS(1)=1 P4F IB=1 P4F DO 300 I=1,M P4F HSPIKE(I)=O P4F 300 CONTINUE P4F C P4F C MAIN LOOP - BLOCKS SETUP P4F C P4F DO 350 I=IS,M P4F TP=BLKS(I) P4F LEN=MPl-IP P4F IF (I.LT.M) LEN=BLKS(I+1)-IP P4F IF (LEN.GT.1) GOTO 310 F4F LENTR=LENTR+1 P4F GOTO 350 P4F 310 IF (LENTR.EQ.0) GOTO 320 P4F C P4F C TRIANGULAR BLOCK P4F C P4F IB=IS+1 P?F BLKS(IB)=IP P4F LEN-R=0 P4F C P4F C PROCESS A BUMP P4F C P4F 320 IE=I?+LEN-1 P4F DO 325 J=1,M P4F 325 RC(J)=0 P4F CALL PARTEL (IP,IE,RI,CS,CNTO,ROTN,CE,RC,NZ) P4F IF (LEN.EQ.2) GOTO 330 P4F CALL BUMP (IP,IE,M,NZ,RI,CS,CE,RC,CNTO,RNTO,ROTN, P4F 1 HSPIKE,NSPIKE,IW,NW) P4F GOTO 340 P4F C P4F C SET PARAMETERS FOR BUMP OF SIZE 2 P4F C ?4F 330 NZ=4 P4F NSPIK7=1 P4F J=CNTO(IE) ?4F HSPIKE(J)=IP ?4F C P4F 310 J=CNTO(IP) ?F '-PIK7(J)=NSPIKE ?4F IF (I.EQ.M) GOTO 360 IB=IS+1 BLKS( 73)=IE+1 4 350 CONTINUE ?4F C P4F 360 BLKS(IB+1)= MP ?4F RETURN P4F END P4F A.3 SUBROUTINE PARTEL (JS,JE,RI,CS,CNTO,ROTN,CE,RC,NZ) PART C PART C PARTBL PARTITIONS THE ROW INDICES LIST FOR A SPECIFIED RANGE PART C OF COLUMNS. FOR EACH COLUMN THE LIST IS PARTITIONED SO THAT PART C ALL THE ELEMENTS IN ROWS WHICH ARE ABOVE THE PARTITION ROW PART C -JE- APPEAR FIRST. USUALLY THE SUBROUTINE IS USED TO PARTITION PART C THE COLUMNS OF A BLOCK, SO THAT ALL THE ELEMENTS IN THE BLOCK PART C BELONG TO THE FIRST PART. PART C PART INTEGER Rl(1),CS(1),CE(1),CNTO(1),ROTN(1),RC(1) PART C PART C INPUT ARGUMENTS PART C JS THE POSITION OF THE FIRST COLUMN TO BE PARTITIONED PART C JE THE POSITION OF THE LAST COLUMN, WHICH IS ALSO THE PART C PARTITION ROW. PART C RI( ) ROW INDICES LIST PART C CS( ) POINTERS TO BEGINNING OF COLUMNS PART C CNTO( ) INVERSE PERMUTATION OF THE COLUMNS PART C ROTN( ) PERMUTATION OF THE ROWS PART C PART C OUTPUT ARGUMENTS PART C CE( ) POiNTERS TO THE END OF THE FIRST PART PART C RC( ) FOR EACH ROW IR IN THE FIRST PART, RC(IR) IS THE NUMBER PART C OF NONZEROS WITHIN THE PARTITIONED COLUMNS. PART C NZ TOTAL NUMBER OF NONZEROS IN THE FIRST PART PART C PART C REMARKS PART C 1. THE ARRAY CNTO( ) IS INDEXED BY THE PRESENT POSITIONS IN THE PART C MATRIX OF THE COLUMNS TO BE PARTITIONED. PARI C 2. ALL OTHER ARRAYS ARE INDEXED BY THE ORIGINAL INDICES. FOR PART C EXAMPLE, THE ROW IR IS IN POSITION NR=ROTN(IR), AND IT PART C BELONGS TO THE FIRST PART IF NR.LE.JE, IN WHICH CASE RC(IR) PART C CONTAINS ITS NONZEROS COUNT. PART C PART NZ=0 PART DO 100 J=JS,JE PART IC=CNTO(J) PART K=CS(IC) PART KE=CS(IC+1)-1 ?ART 10 IF(K.GT.KE) GOTO 30 PART 1R=RI(K) PART NR=ROTN(IR) PART TF(NR.GT.JE) GOTO 20 PART RC(IR)=RC(IR)+1 ?ART NZ=NZ+1 PART K=K-.1 ?ART GOTO 10 PART 20 RI(K)=RI(KE) ?ART RI(KE)=IR PART KE=KE-1 PART GOTO 10 ?ART 30 CS(IC)=KE PART 100 CONTINUE PART RETURN PART END PART A. 4 SUBROUTINE BUMP (JS,JE,M,NZ,RI,CS,CE,RC,CNTO,RNTO,ROTN, BUMP 1 HSPIKE,NSPIKE,IN,NW) BUMP C BUMP C BUMP IS A LINKING SUBROUTINE WHICH ASSIGNS THE BUM? BUMP C PROCESSING TO ONE OF TWO SUBROUTINES: HRF OR HRS. BUMP C BOTH PERFORM THE HR ROUTINE ON A SINGLE BUMP. HRF IS BUMP C FASTER, BUT IT REQUIRES MORE STORAGE SPACE THAN HRS. BUMP C A WORK ARRAY IW( ) OF LENGTH NW IS SUPPLIED. IF THE BUMP C WORK ARRAYS REQUIRED BY HRF CAN FIT INTO IW{ ), THEN BUMP C HRF IS USED. OTHERWISE HRS IS CALLED. BOTH SUBROUTINES BUMP C ASSUME THAT THE BUMP IS IRREDUCIBLE, AND THEY WILL NOT BUMP C IDENTIFY A TRIANGULAR BLOCK. BUMP C BUMP INTEGER RI(1),CS(1),CE(1),CNTO(1),RNTO(1),ROTN(1),HSPIKE(1) BUMP INTEGER IW(NW) BUMP C BUMP C INPUT ARGUMENTS BUMP/A C IW WORK ARRAY OF SIZE NW. USED BY SUBROUTINE HRF IF BUMP/A C NW .GT. 4M + NZ + 2*(JE-JS+2) . BUMP/A C IN THIS CASE THE FIRST M ENTRIES OF IW( ) MUST CONTAIN BURP/A C THE COLUMN PERMUTATION, WHICH IS THE INVERSE OF CNTO( ). BUMP/A C FOR DESCRIPTION OF OTHER ARGUMENTS REFER TO SUBROUTINE HRS. BUMP C BUM? C SUBROUTINES REQUIRED BUMP C ROWREP TO CREATE A ROW-WISE REPRESENTATION OF THE BUMP. BUMP C (USED ONLY WHEN HRF IS CALLED) FuXP C HRF THE FAST HELLERMAN-RARICK ROUTINE BUMP HRS THE HELLERMAN-RARICK ROUTINE WITH MINIMUM STORAGE BUMP C REQUIREMENTS. BUMP C BUMP C BUMP L=JE-JS+2 BUMF IF (4*M.?LL+NZ.GE.NW) GOTO 10 BUMP IA=1+M BUMP IB=IA+NZ BUMP CALL ROWREP (JS,JE,M,RI,CS,CE,RC,CNTO,IW(IA),IW(IB)) BUMP IC=IB+M+1 BUMP ID=IC+M BUMP IE=ID+L BUMP IF=IE+M BUMP CALL HRF (JS,JE,RI,CS,CE,RC,CNTO,IW,RNTO,ROTN,IW(IA),IW(IB), BUMP 1 HSPIKE,NSPIKE,IW(IC),IW(ID),IW(IE),IW(IF)) BUMP RETURN BUMP 10 CALL HRS (JS,JE,RI,CS,CE,RC,CNTO,RNTO,ROTN,HSPIKE,NSPIKE) BUMP RETURN BUMP END BUMP A.5 SUBROUTINE ROWREP (JS,JE,M,RI,CS,CE,RC,CNTO,CI,RS) ROWREP C ROWREP C GIVEN A PARTITIONED COLUMN-WISE REPRESENTATION OF A SPARSE ROWREP C MATRIX AND A PERMUTATION OF ITS COLUMNS, ROWREP CREATES THE ROWREP C ROW-WISE REPRESENTATION OF A PRINCIPAL SUBMATRIX OF THE ROWREP C PERMUTED KATRIX. THIS SUBMATRIX IS COMPOSED OF A SPECIFIED ROhvREP C RANGE OF COLUMNS ( WITHIN THE NEW ORDER ), AND THE ROWS IN ROWREP C THE UPPER PART OF THE PARTITIONED COLUMNS. SINCE BOTH ROWS ROWREP C AND COLUMNS MAY BE PERMUTED, THE ROWS AND COLUMNS OF THE ROWREP C SUBMATRIX DO NOT HAVE THE SAME INDICES. ROWREP C ROWREP INTEGER RI(1),CS(1),CE(1),RC(1),CNTO(1),CI(1),RS(2) ROWREP C ROWREP C INPUT ARGUMENTS ROWREP C JS POSITION OF THE FIRST COLUMN IN THE SUBMATRI. ROWREP C JE POSITION OF THE LAST COLUMPI IN THE SUBMATRIX ROWREP C M ORDER OF THE MATRIX ROWREP C FOR DESCRIPTION OF OTHER ARGUMENTS REFER TO SUBROUTINE HRS ROWREP C ROWREP C OUTPUT ARGUMENTS ROWREP C CI( ) THE ROW-WISE COLUMN INDICES LIST WITHIN THE SUBMATRIX. ROWREP C RS( ) POINTERS TO THE BEGINNING OF ROWS IN CI( ). ROWREP C RS( ) HAS LENGTH M+1, AND IT REFERS TO A ROW BY ITS ROWREP C ORIGINAL INDEX. SO IF ROW I IS NOT IN THE SUBMATRIX ROWREP C RS(I)=RS(1+1). ROWREP C ROWREP C ROWREP M1=M+1 ROWREP RS(1)=1 ROWREP RS(2)=1 ROWREP DO 10 I=3,M1 ROWREP 10 RS(I)=RS(I-1)+RC(I-2) ROWREP DO 30 J=JS,JE ROWREP IC=CNTO(J) ROWREP KS=CS(IC) ROWREP KE=CE(IC) ROWREP DO 20 K=KS,KE ROWREP I=RI(K)+1 ROWREP KK=RS(I) ROWREP CI(KK)=IC ROWREP 20 RS(I)=RS(I)+1 RO'4REP 30 CONTINUE ROWREP RETURN ROWREP END ROWREP A.6 SUBROUTINE HRF (JS,JE,RI,CS,CE,RC,CNTO,COTN,RNTO,ROTN,CI,RS, HRF 1 HSPIKE,NSPIKE,CC,TH,TF,LH) HRF C zRF C HRF IS THE HELLERMAN-RARICK ROUTINE ADJUSTED TO OPERATE ON HRF C AN IRREDUCIBLE BLOCK (BUMP). IT USES A PARTITIONED COLUMN- HRF C WISE REPRESENTATION AS WELL AS A LOCAL ROW-WISE ARRAYS. FREE HRF C COLUMNS AND R')WS ARE KEPT SORTED IN ORDER TO REDUCE EXECUTION HRF C TIME. THE SUBROUTINE EXPECTS A SQUARE IRREDUCIBLE MATRIX, HRF C AND IT REARRANGES THE ROWS AND COLUMNS SO AS TO MINIMIZE THE HRF C NUMBER OF SPIKES. HRF C HRF INTEGER RI(1),CS(1),CE(1),COTN(1),CNTO(1),ROTN(1),RNTO(1),CI(1) HRF INTEGER RS(1),CC(1),RC(1),TH(1),TF(1),LH(1),HSPIKE(1) ERF C KRF C INPUT ARGUMENTS HRF C JS POSITION OF THE FIRST COLUMN IN THE BUMP :;RF C JE POSITION OF THE LAST COLUMN IN THE BUMP HRF C RI( ) COLUMN-WISE ROW INDICES LIST ;iRF C CS( ) POINTERS TO THE BEGINNING OF COLUMN IN RI( ) HRF C CE( ) POINTERS TO END OF THE BUMP IN RI( ) HRF C RC( ) ROWS NONZEROS COUNTS WITHIN THE BUMP ;;R7 C CNTO( ) INITIAL INVERSE COLUMN PERMUTATION. CNTO(I) IS THE HRF C ORIGINAL INDEX OF THE COLUMN IN THE I'TH POSITION. HF? C COTN( ) INITIAL COLUMN PERMUTATION HRF C RNTO( ) INITIAL INVERSE ROW PERMUTATION F_RF C ROTN( ) INITIAL ROW PERMUTATION HRF C CI( ) ROW-WISE COLUMN INDICES LIST WITH:N THE BUMP ;RF C RS( ) POINTERS TO BEGINNING OF ROWS IN CI( ERF C HRF C OUTPUT ARGUMENTS HRF C CNTO/RNTO NEW INVERSE COLUMN/ROW PERMUTATION HRF C COTN/ROTN NEW COLUMN/ROW PERMUTATION HRF C HSPIKE( ) IF HSPIKE(J)=0, THEN COLUMN J IS NOT A SPIKE. ER C OTHERWISE, HSPIKE(J) IS THE HEIGHT OF THE SPIKE HRF C NSPIKE NUMBER OF SPILES IN THE BUMP HRF C HRF C INTERNAL USE OF ARGUMENTS HRF C CC( ) COLUMNS NONZEROS COUNTS WITHIN TEE BUMP ;RF C TH( ) COLUMNS ARE ORDERED IN CNTO( ) BY INCREASING TALLY ERF C VALUE. TH(I) POINTS TO THE FIRST COLUMN WITH VALUE 1. FRF C TF( ) TALLY FUNCTION VALUES ERF C LH( ) ROWS ARE ORDERED IN RNTO( ) BY INCREASING NONZEROS ERF C COUNTS. LH(I) POINTS TO THE FIRST ROW WITH COUNT I. ERF C ERF C SUBROUTINES REQUIRED ER? C LSORT TO SORT THE ROWS BY NONZEROS COUNT HRF C MAXTLY TO FIND THE COLUMN WITH MAXIMUM WEIGHTED TALLY VALUE HRF C REDRC TO REDUCE ROW COUNTS WHENEVER A COLUMN IS ASSIGNED C TALLY TO COMPUTE OR UPDATE TALLY FUNCTION VALUES HF C ERF * **** C ERF C SET INITIAL VALUES RF C HR? NU=JS HRF L=JE hRF NSPIKE=0 IRF A.7 C HRF 59 DO 5 J=JS,jE hRF 60 IC=CNTO(J) HRF 61 5 CC(IC)=CE(IC)-CS(IC)+1 HRF 62 C HRF 63 C FORM ROW LIST SORTED BY ROW COUNTS HRF 61, C HRF 65 CALL LSORT (JS,JE,RNTO,ROTN,RC,LH,TF) HRF 66 C HRF 67 C SET VALUES FOR INITIAL SPIKE SELECTION HRF 68 C HRF 69 TH(1)=NU HRF 70 IR=RNTO(NU) HRF 71 MINRC=RC(IR) HRF 72 KE=LH(MINRC+1)-1 HRF 73 C HRF 74 C COMPUTE TALLY FUNCTION FOR ALL FREE COLUMNS HRF 75 C HRF 76 10 KL=KE-NU+2 HRF 77 IS=TH(1) HRF 78 DO 20 I=IS,L HRF 79 IC=CNTO(I) HRF 80 TF(IC)=O HRF 81 20 CONTINUE HRF 82 C HRF 83 IPCH=L+1 HRF 84 DO 30 I=1,KL HRF 85 TH(I)=IPCH HRF 86 30 CONTINUE HRF 87 C HRF 88 CALL TALLY (NU,KE,1.CI,RS,CC,TF,TH,CNTO,COT,RNTO) HRF 89 C HRF 90 C FIND HOW MANY COLUMNS WITH MA..MUM TALLY HRF 91 C HRF 92 40 ICH=CNTO(L) HRF 93 ITALLY=TF(ICH) HRF 94 IS=TH(ITALLY) HRF 9c IP=L HRF 9E IC=ICH HRF 97 IF (IS.EQ.L.OR.KE.EQ.JE) GOTO 80 HRF 98 C HRF 9; C SELECT ONE FROM THE MAXIMUM TALLY COLUMNS... HRF 10C C HRF 101 IF (ITALLY.GT.1) GOTO 50 HRF 102 C HRF 10 C ...BY WEIGHTED TALLY (IF MAXIMUM TALLY = 1) HRF 10L C HRF 10- CALL MAXTLY (IS,L,RI,CS,C,RC,CNTO,MINRC,I?) HRF 10 GOTO 70 HRF 10, C HRF 10, C ... BY MAXIMUM COLUMN COUNT HRF 1o C HRF lic 50 MAXCC=0 HRF 11 DO 60 I=IS,L HRF 112 IC=CNTO(I) HR 11 IF (MAXCC .GE. CC(IC)) GOTO 60 HRF 1 MAXCC=CC(IC) HRF 11 A.8 IP=I HRF 116 60 CONTINUE HRF 117 70 IC=CNTO(IP) HRF 118 C HRF 119 C CHECK IF A PIVOT CAN BE ASSIGNED HRF 120 C HRF 121 80 IF (MINRC.EQ.1) GOTO 110 HRF 122 C HRF 123 C STACK A SPIKE HRF 124 C HRF 125 IF (IP.EQ.L) GOTO 90 HRF 126 COTN(IC)=L HRF 127 COTN(ICH)=IP HRF 128 CNTO(L)=IC HRF 129 CNTO(IP)=ICH HRF 130 90 CC(IC)=0 HRF 131 L=L-1 HRF 132 HSPIKE(IC)=NU HRF 133 NSPIKE=NSPIKE+1 HRF 134 C HRF 135 C REDUCE ROW COUNTS FOR SPIKE HRF 136 C HRF 137 CALL REDRC (CS(IC),CE(IC),RI,RC,RNTO,ROTN,LH) HRF 133 C HRF 139 C CHECK WHETHER TO RECOMPUTE OR UPDATE TALLY HRF 140 C HRF 141 KEO=KE HRF 142 KE=LH(MINRC)-1 HRF 143 MINRC=MINRC-1 HRF 144 NUP=KEO-KE HRF 145 IF (NUP/(KE-NU+1) .GE. 1) GOTO 10 HRF 146 C HRF 147 C U?DATw TALLY HRF 148 C HRF 149 IF (NUP.EQ.0) GOTO 40 HRF 150 ITS=ITALLY+1 HRF 151 IPCH=L+1 HRF 152 DO 100 I=ITS,KL HRF 153 TH(I)=IPCH HRF 154 100 CONTINUE HRF 155 CALL TALLY (KE+1,KEO,0,CI,RS,CC,TF,TH,CNTO,COTN,RNTO) HRF 1S6 GOTO 40 HRF 157 C HRF 158 C REDUCE ROW COUNTS FOR PIVOT COLUMN HRF 159 C HRF 160 110 CALL REDRC (CS(TC),CE(IC),RI,RC,RNTO,ROTN,Lii) d?.F 161 C HRF 162 C ASSIGN PIVOT OR SPIKE TO POSITION NU HRF 163 C HRF 164 120 IH=CNTO(NU) HRF 165 COTN(IC)=NU HRF 166 COTN(IH)=IP HRF 167 CNTO(NU)=IC HRF 168 CNTO(IP)=IH. HRF 169 CC(IC)=0 HRF 170 A.9 C HRF 171 C TEST FOR TERMINATION HRF 172 C HRF 173 NU=NU+1 HRF 174 IF (NU.GT.L) RETURN HRF 175 C HRF 76 C SCAN ROW COUNTS FOR NEXT STEP HRF 177 C HRF 178 IR=RNTO(NU) HRF 179 MINRC=RC(IR) HRF 180 IF (MINRC .EQ. 0) GOTO 150 HRF 181 KE=LH(MINRC+1)-1 HRF 182 IF (MINRC.GT.1 .OR. KZE.GT.NU) GOTO 10 HRF 183 C HRF 184 C UNIQUE PIVOT ROW - FIND PIVOT COLUMN ERF 185 C HRF 186 IS=RS(IR) HRF 187 IE=RS(IR+1)-1 HRF 188 DO 130 I=IS,IE HRF 189 IC=CI(I) HRF 190 IF (CC(IC).NE.0) GOTO 140 HRF 191 130 CONTINUE HRF 192 140 IP=COTN(IC) HRF 193 GOTO 110 HRF 194 C HRF 195 C UNSTACK A SPIKE HRF 196 C HRF 197 150 L=L+1 HRF 198 IP=L HRF 199 ICmCNTO(L) HRF 200 GOTO 120 HRF 201 C HRF 202 END HRF 203 A.10 SUBROUTINE HRS (JS,JE,RI,CS,CE,RC,CNTO,RNTO,ROTN,HSPIKE,NSPIKE) -RS 2 C HRS 3 C HRS IS A VERSION OF THE HELLERMAN-RARICK ROUTINE WHICH 4RS 4 C OPERATES ON AN IRREDUCIBLE BLOCK (BUMP). IT USES ONLY A =RS C PARTITIONED COLUMN-WISE REPRESENTATION. THE SUBROUTINE :RS 6 C EXPECTS A SQUARE IRREDUCIBLE MATRIX, AND IT REARRANGES HRS 7 C THE ROWS AND COLUMNS SO AS TO MINIMIZE THE NUMBER OF SPIKES. LRS 8 C ERS 9 INTEGER RI(1),CS(1),CE(1),RC(1),CNTO(1),RNTO(1),HSPIKE(1),ROTN(1) HRS .O C HRS 1. C INPUT ARGUMENTS HRS C JS POSITION OF THE FIRST COLUMN IN THE BUMP HRS 13 C JE POSITION OF THE LAST COLUMN IN THE BUMP HIS 14 C RI( ) COLUMNWISE ROW INDICES LIST HRS 15 C CS( ) POINTERS TO BEGINNING OF COLUMNS IN RI( )HRL .6 C CE( ) POINTERS TO END OF THE BUMP IN RI( ) .RS 17 C RC( ) ROWS NONZEROS COUNT WITHIN THE BUMP HRS TS C CNTO( ) INITIAL INVERSE COLUMN PERMUTATION. CNTO(I) IS THE R.S .9 C ORIGINAL INDEX OF THE COLUMN IN THE I'TH POSITION. HPS 20 C RNTO( ) INITIAL INVERSE ROW PERMUTATION HRS 21 C ROTN( ) THE INITIAL ROW PERMUTATION HRS 22 C ERS 23 C OUTPUT ARGUMENTS HRS 24 C CNTO( ) NEW INVERSE COLUMN PERMUTATION HRS 25 C RNTO( ) NEW INVERSE ROW PERMUTATION HRS 26 C ROTN( ) NEW ROW PERMUTATION RS 27 C HSPIKE( ) IF HSPIKE(J)=O, THEN COLUMN J IS NOT A SPIKE. HRS 28 C OTHERWISE HSPIKE(J) IS THE HEIGHT OF THE SPIKE. ERS 29 C NSPIKE NUMBER OF SPIKES IN THE BUMP RS 30 C HRS 31 C SUBROUTINES REQUIRED HRS 32 C MAXTLY TO FIND THE COLUMN WITH MAXIMUM WEIGHTED TALLY VALUE HRS 33 C FRS 34 C REMARKS HRS Z5 C 1. ALL ARRAYS ARE INDEXED WITH RESPECT TO THE WHOLE MATRIX. HRS 35 C FOR EXAMPLE, CNTO(J) IS THE ORIGINAL INDEX OF THE COLUMN WHICH HRS 37 C IS NOW IN THE J'TH POSITION OF THE MATRIX (NOT THE BUMP). HRS 35 C 2. THE CONTENT CF ALL ARRAYS (EXCEPT RC( ) ) REFERS TO THE BIG -RS 39 C MATRIX, IN WHICH THE BUMP IS A DIAGONAL BLOCK. HRS !0 C HRS 41 C HRS 42 C SET INITIAL VALUES HRS C HRS -4 NU=JS r3S 45 L=JE HRS 45 NSPIKE=0 HS 47 C ?RS C FIND MINIMUM ROW COUNT 'Rs C His 50 10 MINRC=JE HRS 5: DO 20 I=NU,JE .RS I IR=RNTO(I) H RS 5: NELEM=RC(IR) ;-S 5 IF (MINRC.GT.NELEM) MINRC=NELEM RS 55 20 CONTINUE EES A. 11 C HRS 57 C FIND A FREE COLUMN WITH VAXIMUM TALLY HRS 58 c iRS 59 30 CALL MAXTLY(NU,L,RI,CS,CE,RC,CNTO,M.INRC-1,IP) HRS 60 C HRS 61 C CHECK IF A PIVOT CAN BE ASSIGNED HRS 62 C HRS 63 IF (MLNRC.EQ.1) GOTO 100 HRS 64 c HRS 65 C STACK A SPIKE HRS 66 C HRS 67 IC=CNTO(IP) HRS 68 CNTO(IP)=CNTO(L) HRS 69 CNTO(L)=IC HRS 70 HSPIKE(IC)=NU HRS 71 ',=L-1 HRS 72 NSPIKE=NSPIKE+1 HRS 73 C HRS 74 C REDUCE ROW COUNTS HRS 75 C HRS 76 IS=CS(IC) -IRS 77 7=CE(IC) HRS 78 40 I=IS,IE HRS 79 IR=RI(I) HRS 80 UO RC(IR)=RC(IR)-1 HRS 81 MINRC=MINRC-1 HRS 82 GOTO 30 HRS 83 C HRS 84 C FORM A CHAIN OF ROWS hITH COUNT 1 HRS 85 C REDUCE ROW COUNTS FOR OTHER FREE ROWS HRS 86 C HRS 87 100 IC=CNTO(IP) HRS 88 MP=0 HRS 89 IS=CS(IC) HRS 90 IE=CE(IC) HRS 91 DO 120 I=IS,IE HRS 92 IR=RI(I) HRS 93 NELEM=RC(IR) HRS 94 IF (NELEM.EQ.0) GOTO 120 HRS 95 IF (NELEM.GT.1) GOTO 110 HIRS 96 RC(IR)=MP HRS 97 MP=IR HRS 98 GOTO 120 HRS 99 110 RC(IR)=NELEM-1 MRS 100 120 CONTINUE HRS 101 C MRS 102 C ASSIGN PIVOT OR SPIKE TO POSITION NU HRS 103 C HRS 104 130 CNTO(IP)=CNTO(NU) .IRS 105 CNTO(NU)=IC :RS 106 C HRS 107 I=ROTN(MP) HRS 108 IR=RNTO(NU) MRS 109 RNTO(I)=IR HRS 110 RNT0(NU)=MP HRS 111 ROTN(IR)=I HRS 112 ROTN(MP)=NU HRS 113 A.12 C HRS 114 C TEST FOR TERMINATION HRS 115 C HRS 116 NU=NU+1 HRS 117 IF (NU.GT.L) GOTO 200 HRS !18 C HRS 119 C SEE IF A SPIKE CAN BE UNSTACKED HRS 120 C HRS 121 IR=RC(MP) HRS *22 IF (IR.EQ.0) GOTO ' J RC(MP)=0 MP=IR H; C HRS C UNSTACK A SPIKE HRS le, C HRS 128 L=L+1 HRS 129 IP=L HRS 130 IC=CNTO(L) HRS 131 GOTO 130 HRS 132 C HRS 133 200 RETURN HRS 134 END HRS 135 A.13 SUBROUTINE LSORT (IS,IE,CNTO,COTN,CC,LH,WORK) LSORT 2 C LSORT 3 C LSORT SORTS A SPEC:FIED RANGE OUT OF A ?F.RMUTATTON VECTOR LSORT 4 C CNTO( ) ACCORDING TO COUNTS GIVEN BY THE VECTOR CC( ). LSORT 5 C THIS MEANS THAT UPON RETURN THE COUNT OF CNTO(I) IS LESS THAN LSORT 6 C OR EQUAL TO THE COUNT OF CNTO(I+1) FOR ANY I BETWEEN IS AND LSORT 7 C IE-1. ALL COUNTS ARE KNOWN TO BE POSITIVE INTEGERS. LSORT 8 C LSORT 9 C LH(I) POINTS TO THE START OF THE LIST OF ELEMENTS WITH COUNT I. !SORT 10 C COTN( ) IS MAINTAINED AS INVERSE OF CNTO( ). LSORT 11 C WORK( ) IS A WORK ARRAY WHICH SHOULD HAVE THr, SAME SIZE AS LSORT 12 C CNTO( ), COTN( ) AND CC( ). LSORT 13 C LSORT 14 C LSORT 15 INTEGER WORK(1),COTN(1),CNTO(1),CC(),LH(1) LSORT 16 C LSORT 17 C INITIALIZE LISTHEADS TO BE NIL LSORT 18 C LSORT 19 M=IE-IS+2 LSO : 20 DO 10 I=1,M LSORT 21 LH(I)=O LSORT 22 10 CONTINUE LSORT 23 C LSORT 24 C GENERATE LITS OF EQUAL COUNT USING WORK LSORT 25 C LSORT 26 DO 20 I=IS,IE LSORT 27 II=CNTO(I) LSORT 28 L=CC(II) LSORT 29 WORK(II)=LH(L) LSORT 30 LH(L)=II LSORT 31 20 CONTINUE LSORT 32 C LSORT 33 C REPLACE INTO CNTO AND SAVE START OF LISTS LSORT 34 C LSORT 35 I=IS LSORT 36 DO 40 NUM=1,M LSORT 37 IP=LH(NUM) LSORT 38 LH(NUM)=I LSORT 39 30 IF (IP.EQ.0) GOTO 40 LSORT 40 COTN(IP)=I LSORT 41 CNTO(I)=I? LSORT D2 IP=WORK(IP) LSORT 4, I=I+1 LSORT '.4 GOTO 30 LSORT 45 40 CONTINUE LSORT 46 RETURN LSORT 47 END LSORT 48 A.14 SUBROUTINE TALLY (KS,KE,OPER,CI,RS,CC,TF,TH,CNTO,COTN,RNTO) TALLY 2 C TALLY 3 C TALLY COV-PUTES OR UPDATES TALLY FUNCTION VALUES. TALLY 4 C THE TALLY FUNCTION IS THE NUMBER OF NONZEROS A COLUMN HAS IN TALLY 5 C A SPECIFIED RANGE OF ROWS. USUALLY THIS RANGE INCLUDES ALL TALLY 6 C THE ROWS WITH MINIMAL NONZEROS COUNT. THE COLUMNS ARE KEPT TALLY 7 C IN A LIST SORTED BY INCREASING TALLY VALUES. TALLY 8 C TALLY 9 IPITFGER CI(1),RS(1),CC(1),TF(1),TH(1),COTN(1),CNTO(1),RNTO(1) TALLY 10 INTEGER OPER TALLY 11 C TALLY 12 C INPUT ARGUMENTS TALLY 13 C KS POSITION IN RNTO( ) OF THE FIRST ROW TO BE SCANNED TALLY 14 C KE POSITION OF THE LAST ROW TALLY 15 C OPER IF OPER=1 TALLY VALUES ARE COMPUTED BY ADDING I TALLY 16 C TO TF(IC) WHENEVER A NONZERO IS FOUND IN COLUMN IC. TALLY 17 C IF JPER=O TALLY VALUES ARE UPDATED BY SUBTRACTING TALLY 18 C I FROM TF(IC). TALLY 19 C FOR DI;RlPTION OF OTHER ARGUMENTS REFER TO SUBROUTINE HRF. TALLY 20 c TALLY 21 C TALLY 22 DO 30 KK=KS,KE TALLY 23 II=RNTO(KK) TALLY 24 IS=RS(II) TALLY 25 IE=RS(IIl1)-i TALLY 26 C TALLY 27 DO 20 I=IS,IE TALLY 28 IC=CI(I) TALLY 29 IF (CC(IC) .LE. 0) GOTO 20 TALLY 30 K=TF(IC)+OPER TALLY 31 TF(IC)=K-1+OPER TALLY 32 L=TH(K)-OPER TALLY 33 IP=COTN(IC) TALLY 34 TH(K)=L+1-OFER TALLY 35 IH=CNTO(L) TALLY 36 COTN(IC)=L TALLY 37 COTN(IH)=IP TALLY 38 CNTO(L)=IC TALLY 39 CNTO(IP)=IH TALLY 40 20 CONTINUE TALLY 41 C TALL' 42 30 CONTINUE TALLY 43 RETURN TALLY 44 END TALLY 45 A.15 SUBROUTINE MAXTLY(JS,JE,R1,CS,CE,RC,CNTO,IOFF,JP) MAXTLY 2 C MAXTLY 3 C COMPUTE MAXIMUM WEIGHTED TALLY AMONG ALL RELEVANT COLUMNS MP-XTLY 4 C MAXTLY 5 INTEGER RI(1),CS(1),CE(0),RC(1),CNTO(1) MAXTLY 6 DIMENSION TVAL(10) MAXTLY 7 DATA TVAL/E25,1E17,IE15,1E13,1E11,1E9,1E7,1E5,1E3,1./ MAXTLY 8 C MAXTLY 9 C INPUT ARGUMENTS MAXTLY 10 C JS THE FIRST COLUMN TO BE EVALUATED MAXTLY 11 C JE THE LAST CO! UMN MAXTLY 12 C RC( ) THE CURRENT COUNT OF REMAINING NONZEROS MAXTLY 13 C IOFF. A CONSTANT TO BE SUBTRACTED FROM THE NONZEROS COUNT. MAXTLY 14 C MAY BE USED FOR EXCLUDING ROWS WITH LOW COUNT. MAXTLY 15 C FOR DESCRIPTION OF OTHER INPUT ARGUMENTS SEE SPIKE. 1-iAXTLY 16 C MAXTLY 17 C OUTPUT ARGUMENTS MAXTLY 18 C JP THE COLUMN WITH MAXIMUM TALLY MAXTLY 19 C MAXTLY 20 TMAX=-1. MAXTLY 21 DO 20 J=JS,JE MAXTLY 22 IC=CNTO(J) MAXTLY 23 T=0. MAXTLY 24 KS=CS(IC) MAXTLY 25 KE=CE(IC) MAXTLY 26 DO 10 K=KS,KE MAXTLY 27 IR=RI(K) MAXTLY 28 NZ=RC(IR)-IOFF MAXTLY 29 IF(NZ.GT.10) NZ=10 MAXTLY 30 IF(NZ.GT.0 ) T=T+TVAL(NZ) MAXTLY 31 10 CONTINUE MAXTLY 32 IF(TMAX.GE.T) GOTO 20 MAXTLY 33 TMAX=T MAXTLY 34 JP=J MAXTLY 35 20 CONTINUE MAXTLY 36 RETURN MAXTLY 37 END MAXTLY 38 A.16 SUBROUTINE REDRC (IS,IE,RI,RC,RNTO,ROTN,LP) REDRC 2 C REDRC 3 C REDRC REDUCES RO COUNTS AND MAINTAINS A SORTED ROWS LIST. REDRC 4 c REDRC 5 INTEGER fI(1),RC(1),ROTN(1),RNTO(1),LH(1) REDRC 6 C REDRC 7 C INPUT ARGUMENTS REDRC 8 C IS CTART POINTER TO RI( ) REDRC 9 C IE END POINTER TO RI( ). THE ROWS WHOSE COUNT HAS TO BE REDRC 10 C REDUCED ARE ON THE ROW INDICES LIST FROM RI(IS) TO RI(IE) REDRC 11 C RI( ) THE ROW INDICES LIST REDRC 12 C REDRC 13 C INPUT/OUTPUT ARGUMENTS REDRC 14 C RC( ) ROWS NONZEROS COUNTS REDRC 15 C RNTO( ) THE INVERSE ROW PERMUTATION WHICH CONTAINS THE REDRC 16 C SORTED ROWS LIST REDRC 17 C ( ) THE ROW PERMUTATION REDRC 18 C Lt. ,) POINTERS TO THE ROWS LIST. LH(I) POiNTS TO THE REDRC 19 C FIRST ROW WITH NONZERO3 COUNT I. REDRC 20 C REDRC 21 C REDRC 22 DO 10 l=IS,IE REDRC 23 IR=RI(I) REDRC 24 K=RC(IR) REDRC 25 IF (K .LE. 0) GOTO 10 REDRC 26 RC(IR)=K-1 REDRC 27 L=LH(K) REDRC 28 IP=ROTN(IR) REDRC 29 LH(K)=L+1 REDRC 30 IH=RNTO(L) REDRC 31 ROTN(IR)=L REDRC 32 ROTN(IH)=IP REDRC 33 RNTO(L)=IR REDRC 34 RNTO(IP)=IH REDRC 35 10 CONTINUE REDRC 36 RETURN REDRC 37 END REDRC 38