RANMAT - A Random Matrix Generator with Structure Control by Yonatan Levy and Alexander Meeraus Technical Note No. 9 Research Project 671-58 November 1978 Development Research Center World Bank t 1818 H Street, N.W. Washington, D.C. 20433 Preliminary and Confidential: Not for quotation or attribution without prior clearance from the authors. Views expressed are those of the authors and do not necessarily reflect those of the World Bank. Intrc'uction Anyone with the experience of testing programs dealing with matrices or graphs, has felt the need for a matrix generator, and espec4ally one that can generate a variety of matrices. Yet there are no published algorithms (as far as we know) which can carry out this task. This is surprising, consider- ing the many algorithms for random generation of distributions and permuta- tions which are available. Our work offers a random matrix generator which is aimed particularly at sparse matrix generation. RANMAT is a Fortran subroutine which enables the user to specify the sparseness as well as the block lower triangular structure of a square matrix. The structure, however, is hidden by permuting the rows. Subject to the user's specifications, RAN,-AT randomly assigns locations for the nonzero elements, and it may generate numerical values too. The description wLich follows has four sections. The first one describes the options and the input parameters, and it should suffice anyone who just wants to use the subroutine. Section two discusses the principles of the algorithm, and section three deals with implementation details.. Section four contains an illustration and performance details. 1. Ootions and Control Parameters The program offers three options concerning the numerical values of te matrix elements. They may not be generated at all (IV = 0), and in this case only the locations of the nonzero elements are determined. This option is particularly suitable for graphs. The other two options generate val-es for the nonzeros. With one (IV = 1) all the values are between .1 and 1, while the other option (IV = 2) generates an invertible - 2 - matrix by making the (permuted) diagonal elements dominant. In any option the generated matrix is structurally nonsingular. which means that a nonzero diagonal can be obtained by permuting either the rows or the columns. The generated matrix is actually a transformation (by a permutation matrix) of an underlying matrix which is Block Lower Triangular (BLT). The algorithm operates on the latter matrix, and the user can control the sparseness and the BLT form through four parameters: AVER, STD, PERTR and NB. We now briefly describe the role of each one of them. The sparseness is determined by AVER, the averpge number of nonzeros per colLmn. This means that the total number of nonzeros will be approximately M*AVER, where M is the order of the matrix. The exact number is returned by N. The BLT form implies that there are two types of columns: The columns which are contained in blocks, and columns that lie in triangular parts between the blocks. Only "block columns" are allowed to have nonzeros above the diagonal. The requested number of blocks is NB, and the desired percentage of columns in triangular parts is PERT. By varying the values of these two parameters, different. structure patterns can be .obtained, as illustrated in figure 1. We point out that the blocks are reducible, and the number of irreducible blocks may be higher or lower than N3. in the sequel we refer to a block or a triangular segment as a group of columns.- Finally, the variation in the number of nonzeros between columns in the same group, and alse in the size between groups is governed by STD. It is usec' ..fter a proper scaling, as standard deviation for the zCrmal distribution. To summarize, AVER is the average number of nonzeros per column, STD is the standard deviation for all size determinations, PERTR. is the -3- (a) NB =1 (b) NB =6 (c) NB =4 PERTR 100 PERTR =50 PERTR =0 II (d) NB 2 (e) .B 1) NB = PERTR = 30 PERTR = 20 PERTR = 0 F4igure 1: The control of the under!-yin 'g struc::..re by the para:eters \B a3d PERTR. (Nonzeros are assigned to the shaded area only) percentage of "triangular columns", and NB is the number of blocks. The user should note that the actual measures of the generated matrix and the input parameters may somewhat differ. The matrix is written on the file IOUT. The first record contains the size of the matrix, and then each record contains :he row index, column index, and possibly the value of a nonzero element. 2. Princiies of the Algorithm The basic idea of the algorithm is to generate a block lower triangu;'r matrix, column by column (or an upper triangular matrix, row by row), and then to permute the rows (columns) in order to cover the structure. co make it easier for the reader, 'and to avoid repetions, we use, in the following descrption, the variable names from the 7rogra=. The ?rogram has a short preparatory stage, where the mean and standard deviation (s.d.) of the group sizes are determined, and the row ;ermutation is generated. ?TR=PERTR/100 is the portion of columns in :riangular parts. So the average size of a triangular segment is set to D = M*PTR/NB. L'e recall that 1 is the total number of :oLuns, and N3 Is the number of blocks. Similarly the average number of coiuns in a block is 3 7M("1-?TR)!N3. The corresr-nding standard devia:ions are SDD = S=D*D/.VER for triangular parts, and SD3 = S7D*B/AVER for blocks. Note that if STD = 0, -l the blocks are of size 3, and :here are exactl- X3 b'ocks. The main -art of the algorithm generates :e -a:rix as it advances r0o coIumn to column and si-ultaneously from groun : grup. -.s it was -xDlained in the :revious section the columns are groupec, and each column *elon2s either to a block or to a criangular nart. -rouglrout the program I - 5 - is the column currently beins nrocessed. ITR indicates whether the present group is a block (ITR - 0) c- : riangular segment (ITR = 1), and IS is the starting position of the next group. NB is the number of blocks created so far, and N counts the nonzeros. (Initially, I = IS = 1, and ITR - NB = N = 0). Once a column is the first in a group (I = IS), the parameters for this group are calculated as follows. 1. The size of the group, LEN, is evaluated as a normal observation using the mean and s.d. B and SDB for blocks, or D and SDD for triangular parts. 2. To.e average number of nonzeros for columns in the groups, IM. This average is calculated by the formula KD - AVER*(2 - IS +IE - I), where IS and IE are the first and last columns in the group. Note that the mean is adjusted according to the location of the group. Also, if the whole matrix is one group, then IS = 1, IE M, so 3M = AVER, as it should be. 3. The corresponding s.d. is S = STD*)i/AVER. 4. For a block, the height MM of all its columns is the same, and it is equal to the height of the first diagonal element of the block. For example, if a block includes the columns 5 through 10, then LEN = 6, and = M - 4, which means that elements in these columns may be assigned from row 5 and lower. The height of a triangular column I Jq M + I - I. Next, every column in the group is generated in the following manner. 1. The number of nonzeros, K, is det6rmined using the normal distribution with mean XX and s.d. 3. 2. The diagonal element is assigned, and its value is randomly generated if IV = 1, or it is set to K if dialnal do-tinance has been requested. - 6 - 3. The rest K-1 elements are randomly located among the remainder MM-1 positions. Random values In (.1, 1) are assigned if IV is positive. The algorithm terminates after column M has been generated. 3. Implementation Details The main subroutine RANMAT calls three small subprograms: The subroutine SCRAM, which generates a random permutation, and the functions 1UNIF and RNORML, which generate the integer uniform and the normal distribu- :ions respectively. All the programs do not: use COMMON, and all the informa- tion is passed through the argument list. The programs are written in Fortran, and they comply to the ANSI standards except for the use of mixed type expressions. We also use the random number function RANF and the initialization subroutine RANSET, which are CDC utility subprograms. They may have to be replaced for usage on other =achines. 4. Performance The authors have used R.NMAT extensively for testing programs dealing with sp; matrices. Since the set-up work in this program is minimal, the execution time depends linearly,on the number of nonzeros. On the CDC Cyber 73, RANHAT generates abcut 550 elements per second, including Output time. We conclude with an ill-unstration (figure 2). Picture (a) shows a matrix as it was generated by ?.ANYAT. The underlying structure of this matrix is shown in picture (b) with per=tuted row numbers. The third picture is an irreducible BLT form of the sa=e =atrix obtained by both row and column per- mutations. 1111111l11� 111111i11l� 1 11f111111� 12�а5ь7Н901г3ц5b7ду0 I21цSо7890123Ч5ь7q90 1?3u5Ab7143?4деи57у0 1 х,хх...х.х...Х..х;,х ц х;х;� .............•• ц х.х.1.г.............. 2 х х , , , , 18 хх , , , . 18 хх ! . . • г 3 . х . х , . х х. 5 хххх , . , . 5 хх�Х� , . , . � к х , , � � 2 х х, , . . 2 х г х; г , , . S хххх ' 20 х хХ, , . . 20 х хх, , . . b х•,хххкххх...,х.х;г. !2 х..х.х.Х�............ 1 .х...х...,.......... 7•х • х . . . 15 хххх ххх( , . г 12 х х хХ . . • д , х хх х , � � , 7 ,х , xI , , , 15 ххХх ххх , • . Ч х х , хкх -, , 13 . х хххх , г 1Т х х, хх , , . 10 ,� хх хххкх х . 9 х х , ххх , , 13 1 ,х xx,xxj , , 11 ..г....хх...х.....+г 17 х...х.хг..х..,...... 11 .....х.г.дх..�....;.. 12 х х к,с , , , iu , ххх, х,хх , г 1ц t ххкк хх ,, . !s , х хххх , , ti . , хх , х� . . 4 х х , х .хх; , , г iЧ � ххх, х,хх � • 1А , ххх х х,ххх х . В, х,хх х, к, , � 15 хххх ххх , ,. , 1ь , хк , Х, , 10 .х х,х х ххх хх'+ . � !h ....Хх........х..... 10 .х..Хх.,.ххххх,х',.г 19 ...хххх...ххх,ххl.;.. 17 х х,х х • , Ь х ххххххх, х,х . iь , х,х . ,х . 1и хх , , , ,� в, х хх х, , х г ь х ххххх х. х,хх , �g . ххх х х,ххх х , j х х х х 3 х х х . . . . . . . . . х г г0 х„ хх ............:.. 1 х.хх.г.х.х...х..х;.х ! х,хх.х......х..х.х.к Мп г0 �17= ЧТ (�) 'Сlге гпцСr[к дs [t la �eneraCeil. мЕАN= 5.0 SТП= 1.0 (с) '1'he lrrec�uclLle ?,1�ck trtangular form. t�N= 5 PERTp= 1 0• 0 (3 lrreducible blocks) (i,) 'Che underlyin� matrlx, and [tы 1�пrдmetera (4 blockн generated) 1�1��иге 2 Appendix A LISTING OF RAINMAT AND ITS SUBPROGRAMS A. 1 SUBROUTINE RANMAT(M,AVER,STD,PERTR,NB,IV,IST,ICUT,RNTO,N,IW) RANMAT 2 C RANMAT 3 C RANMAT RANDOMLY ASSIGNS LOCATIONS, AND OPTIONALLY GENERATES RANMAT 4 C VALUES, FOR THE NONZERO ELEMENTS OF A SQUARE MATRIX. RANMAT 5 C THE USER SPECIFIES THE DESIRED SPARSENESS AND STRUCTURE OF RANMAT 6 C THE MATRIX. THE PROGRAM GENERATES A BLOCK LCIRER TRIANGULAR RANMAT 7 C MATRIX WITH NONZERO DIAGONAL, AND THE ROWS ARE PERMUTED SO RANMAT 8 C THE STRUCTURE IS NOT RECOGNIZED BY INSPECTION. RANM4T 9 C THE CONSTRUCTION IS COLUMNWISE: FOR EACH COLUMN, ITS HEIGHT RANMAT 10 C AND NUMBER OF NONZEROS ARE DETERMINED, THEN THE NONZEROS ARE RANMAT 11 C RANDOMLY LOCATED, AND VALUES ARE GENERATED IF REQUESTED. RANMAT 12 C RANMAT 13 C OUTPUT RANMAT 14 C THE MATRIX IS WRITTEN COLUMNWISE ON UNIT IOUT. THE FIRST RECORD RANMAT 15 C CONTAINS M, THE ORDER OF THE MATRIX. THEN EACH RECORD CONTAINS RANMAT 16 C THE ROW INDEX, COLUMN INDEX, AND THE VALUE, IF IV IS POSITIVE. RANMAT 17 C (FORMAT IS 215,F15.9). RANMAT 18 C RANMAT 19 INTEGER RNTO(1),1W(1) RkNMAT 20 C RANMAT 21 C VARIABLES (FOR ARRAYS, THE SUBSCRIPTS ARE THE REQUIRED SIZE). RANMAT 22 C RANMAT 23 C INPUT ARGUMENTS RANMAT 24 C M ORDER CF THE MATRIX. (M.LT.100000 OR CHANGE FORMAT) RANMAT 25 C AVER AVERAGE NUMBER OF NONZEROS ?ER COLUMN. THIS AVERAGE IS RANMAT 26 C WEIGHTED ACCORDING TO THE CCLUMN LOCATION. RANMAT 27 C STD STANDARD DEVIATION. USED CO AMOUNT OF NONZEROS PER RANMAT 28 C COLUMN, AND FOR BLOCKS LENGTH. RANMAT 29 C PERTR DESIRED PERCENTAGE OF DIAGONAL ELEMENTS NOT IN BLOCKS. RANMAT 30 C NB REQUESTED NUMBER OF BLOCKS. THE ACTUAL NUMBER MAY VARY. RANMAT 31 C IV IF IV=0 NO VALUES ARE GENERATED, RANMAT 32 C IF IV=1 REAL RANDGM VALUES IN (.1,1) TO ALL NONZEROS, RANMAT 33 C IF IV=2 NONSINGULARITY IS ASSURED BY DIAGONAL DOMINANCE. RANMAT 34 C IST STARTING POINT FOR THE RANDOM NUMBERS GENERATOR - RANMAT 35 C RANSET(IST) - MAY BE USED TO REGENERATE MATRICES. RANMAT 16 C IOUT UNIT NUMBER FOR THE OUTPUT FILE. RANMAT 37 C RANMAT 38 C OUTPUT ARGUMENTS RANMAT 39 C RNTO(M) THE ROW PERMUTATION WHICH MAY BE USED TO REVEAL RANMAT 40 THE STRUCTURE. RANMAT 41 C N TOTAL NUMBER OF NONZEROS. RANMAT 42 C N3 NUMBER OF BLOCKS GENERATED. RANMAT 43 C F.NMAT 44 C INTERNAL USE CF ARGUMENTS RANMAT 45 C Id(M) USED TO MARK THE NGNZERCS adE,.Y LOCATEG. ?ANMAT 46 C RANMAT 47 C RANMAT 46 C SUBPROGRAMS REQUIRF) RANMAT 49 C SCRAm GENERATES A RANDOM PERMUTATION (SUBROUTINE). RANMAT 50 C IUNIF INTEGER RANDOM NUMBERS GENERATCR (FUNCTION). RANMAT 51 C RNRML GENEF'TES AN OBSERVATION FROM NORMAL DISTRIBUTION (FUN). RANMAT 52 C RANMAT 53 C RANMAT 54 C REMARKS RANMAT 55 C 1. ThE LARGEST NJMER TO =E STCED IN ALL ARRAYS IS M. RANMAT 56 C 2. TO GET A MATRIX WITHOUT A S?ECIFIC STRUCTURE, ?UT: NE=1 RANMAT 57 C PERTR=O . RANMAT 58 A.2 C TO GET A COMPLETE TRIASGULAR MATRIX, PUT: NB=1, PERTR=100. RANMAT 59 C RANMAT 60 Cmeeu#*te*A4a*AIm*sIAeaa***eAuessmase*semsseussesseeoIIesemIIuuuIuIuI RANMAT 61 c RANMAT 62 C SET UP PERUTATION AND CCNSTANTS RANMAT 63 C RANMAT 64 CALL SCRAM(RNTO,M) RANMAT 65 CALL RANSET(IST) RANMAT 66 MP1=M+1 RANMAT 67 WRITE(IoUr,90) M,M RANMAT 66 c RANMAT 69 C CALCULATE PARAMJETERS FOR BLOCKS AND TRIANGULAR PARTS RANMAT 70 C RANMAT 71 PTR=PERTR/100. RANMAT 72 B=M*(1.-PTR)/N3 RANMAT 73 D=M*PTR/NB RANMAT 74 SDB=.. RANMAT 75 SDD=0. RANMAT 76 IF (NB.EQ.1) G0 10 RANMAT 77 SDB=STD'B/AVER RAhMAT 18 SDD=STDOD/AVER RANMAT 79 C RAIMAT 80 C INITIALIZE RANMAT 81 c RANMAT 82 10 IS=1 RANMAT 83 ITR=0 RANMAT 84 NB=0 RANNAT 85 N=0 RANMAT 86 FANMAT 7 C MAIN CCLUMNS LOOP RANMAT 88 C RANMAT 89 DO 200 1=1,M RANMAT 90 I7 (IS.GT.I) GOTO 50 RANMAT 91 C SET UP A N-i GROUP RANMAT )2 C RANMAT 93 20 ITR=1-ITR RANMAT 94 IF (ITR.EC.O) GOTO 30 RANJMAT 95 IF (D.EQ.0.) GOTO 20 RANMAT 96 L7N=RNCRM.L(D, SDD) RANMAT 97 IF (LEN) 20,20,40 RANMAT 96 32 LEN=RNCRML(3,SDB) RANNAT 9 IF (LRN.LT.2) AEN=2 ANMAT 100 YA=MP1-IS RANMAT 101 NB=NB+1 RANMAT 102 40 INISLEN RAANMAT 103 IF (IN.GE.M) !N=MP1 RANMAT 1'$ XM=AVER*(2.*M-IS-IN+2./M RANMAT 1c5 S=STD*XM/ AVER RANMAT 106 IS=IN RANMAT 101 RANMAT 10d C DE T ERMINE NUM'ENE3? :NZERCS IN C3LUM! ?ANMAT 109 .,AN,AT 110 50 IF (ITR.EQ.1) MM=M-: ZAiMAT Ill K=RNCRML(XM,S).5 RANMAT '2 IF (K.GT.MM) K=MM RANMAT 113 IF (K.LE.0) K=1 RANMAT 114 DO 60 J=1,MM RANMAT 115 A.3 60 IW(J)=O HANMAT 116 C RANMAT 117 C GENERATE THE PiAGONAL ELEMENT RANMAT lid C RANMAT 119 IIM?1-. RA11MAT 120 IF (1V-1) 110,120,130 RANMAT 121 110 WRITE (IOUT,900) RNTO(II),I RANMAT 122 GOTO 150 RANMAT 123 120 VAL=.9*RANF(,!)+.1 RANMAT 124 GOTO 140 RANMAT 125 130 VAL=K RANMAT 126 140 WRITE (IOUT,900) RNTO(II),I,VAL RANMAT 127 150 IF (K.EQ.1) GOTO 200 RANMAT 128 1l(I1)=1 RANMAT 129 C RASMAT 130 C LOCATE THE NONZEROS HANMAT 131 C RANMAT 132 DO 190 J=2,K RANMAT 133 II=IUNIF(MM) RANMAT 134 163 IF (IW(II).EQ.0) GOTO 170 RANMAT 135 11=11 RANMAT 136 IF (I.GT.MM) II=1 PANMAT 137 GOTu 160 RANMAT 138 170 IF (IV.GT.0) GOTO 180 RANMAT 139 WHITE (IOUT,900) RNTC(II),I RANMAT 140 GOTO 190 RAINMAT 141 180 VAL=.9*RANF(0)+.1 RANMAT 142 WRITE (IOUT,900) RSTC(II),I,VAL RANMAT 143 190 :W(II)=1 RANMAT 144 200 N=N.K RANMAT 145 C RANMAT 146 900 FORMAT(2I5,F15.9) RASMAT 147 RETURN RANMAT 148 END RANMAT 149 A.4 SUBROUTINE SCRAM(NTO,M) SCRAM 2 C SCRAM 3 C SCRAi GENERATES A RANDOM PERMUTATION OF 1,2,...,M. SCRAM 4 C THE PERMUTATION VECTOR IS RETURNED IN NTO( ). SCRAM 5 C SCRAM 6 C EXTERNAL EUNCTION REQUIRED SCRAM 7 C IUNIF INTEGER RANDOM NUMBERS GENERATOR. SCRAM 8 C SCRAM 9 C SCRAM 10 INTEGER NTO(1) SCRAM 11 C SCRAM 12 DO 10 T=1,M SCRAM 13 NTO(I)=I SCRAM 14 10 CONTINUE SCRAM 15 C SCRAM 16 DO 20 I=1,M SCRAM 17 II=IUNIF(M) SCRAM 18 KKnNTO(II) SCRAM 19 NTO(Il)=NTO(I) SCRAM 20 NTO(I)=KK SCRAM 21 20 CONTINUE SCRAM 22 RETURN SCRAM 23 END SCRAM 24 FUNCTION IUNIF(M) IUNIF 2 C IUNIF 3 C IUNIF(M) CHOOSES RANDOMLY ONE NUMBER OUT OF 1,2,...,h . IUNIF 4 C IUNIF 5 IItRANF(0)*M+1 IUNIF 5 IF(II .GT. M) II=M IUNIF 7 IUNIF=II IUNIF 8 RETURN IUNIF 9 END IUNIF 10 FUNCTION RNORML(X,S) RNORML 2 C RNORML 3 C GENERATES ONE C3SERVATI. 'ROM A NORMAL DTSTRIBUTION RNGRML 4 C WITH MEAN X AND STANDARL VIATION S. RNORML 5 C RNOR!ML 6 IF (S.EQ.0.) GOTO 20 RNORML 7 SUM=o. RNORML 8 DO 10 I=1,12 RNORML 9 SUM=SUPl+RANF(0) RNORML 10 10 CONTINUE RNORML 11 20 RNORML=S*(SUM-6.)+X RNORML 12 RETURN RNORML 13 END RNORML 14 Appendix B LISTING OF THREE UTILITY ROUTINES FOR SPARSE MATRICES: RD',T - To read a matrix which is stored in a sparse form. COL?O - To create a row-wise representation of a sparce matrix out of its colur=- wise storage (or vice versa). PRSY' - To generage the symbolic picture of a matrix. (See the picutres in fiugure 2) B.1 SUBROUTINE RDMAT(M,NZ,IV,CS,hi,VLIST,CC,RC,IT) RDMAT 2 C RDMAT 3 C PURPOSE RDMAT 4 C TO READ A SPARSE MATRIX OR GRAPH RDMAT 5 C EITHER WITH OR WITHOUT NUMERICAL VALUES. RDMAT 6 C WE ASSUME THAT THE MATRIX IS STORED IN SPARSE FORM, RDMAT 7 C SO THAT THE FIRST RECORD CONTAIN THE SIZE, AND THEN EACH RECORD RDMAT 8 C REFERS TO ONE NONZERO ELEMENT. IT CONTAINS TWO LOCATION INDICES RDMAT 9 C AND POSSIBLY A VALUE. IT IS ASSUMED THAT THE 1ST INDEX IS RDMAT 10 C NESTED WITHIN THE 2ND, AND THAT THE 2ND INDEX IS NON-DECREASING RDMAT 11 C FOR CONVENIENCE WE REFER TO THE 2ND INDEX AS COLUMN, AND TO THE RDMAT 12 C 1ST AS ROW. RDMAT 13 C RDMAT 14 C RDMAT 15 INTEGER CS('),R!(1),CC(1),RC(1) RDMAT 16 DIMENSION VLIST(1) RDMAT 17 C RDMAT 18 C DESCRIPTION OF VARIABLES (FOR ARRAYS, THE SUBSCIFT RDMAT 19 C IS THE REQIRED SIZE) RDMAT 20 C RDMAT 21 C INPUT ARGUMENTS RDMAT 22 C M DIMENSION LIMITATION ON THE SIZE OF MATRIX. RDMAT 23 C NZ DIMENSION LIMITATION ON THE NUMBER OF NONZEROS. RDMAT 24 C IV IF IV=1, NUMERICAL VALUES ARE EXPECTED FOR THE NONZERO RDMAT 25 C ELEMENTS. RDMAT 26 C IT UNIT NUMBER OF THE INPUT FILE. RDMAT 27 C IF IT=O, DEFAULT VALUE IS 5. RDMAT 28 C RDMAT 29 C OUTPUT ARGUMENTS RDMAT 30 C M SIZE OF THE MATRIX. RDMAT 31 C NZ NUMBER OF NONZEL3 ELEMENTS. RDMAT 32 C CS(M+1) COLUMNS STARTING POINTS IN RI( ). RDMAT 33 C RI(NZ) LIST OF ROW INDICES. RDMAT 34 C VLIST(NZ) LIST OF VALUES. GIVEN ONLY IF IV=1. RDMAT 35 C CC(M) NONZEROS COUNT FOR COLUMNS. RDMAT 36 C RC(M) NONZEROS COUNT FOR ROWS. RDMAT 37 C RDMAT 38 C RDMAT 39 IF(IT.EQ.0) IT=5 RDMAT 40 NBIG=NZ+1 RDMAT 41 MBIG=M RDMAT 42 READ(IT,900) M RDMAT 43 900 FORMAT(215,F15.9) RDMAT 44 IF(M.GT.MBIG) GOTO 800 RDMAT 45 C RDMAT 46 M1=M+1 RDMAT 47 CALL MEMSET(O,CC,M) RDMAT 43 CALL MEMSET(O,RC,M) RDMAT 49 CS(1)=1 RDMAT 50 I=1 RDMAT 51 C RDMAT 52 DO 150 NZ=1,NBIG RDMAT 53 IF(IV.EQ.1) GOTO 100 RDMAT 54 READ(IT,900) IR,IC RDMAT 55 GOTO 105 RDMAT 56 100 READ(IT,900) IR,IC,VLIST(NZ) RDMAT 57 105 1I(EOF(IT)) 180,110 RDMAT 58 B.2 110 CC(IC)=CC(IC)+1 RDMAT 59 RC(IR)=RC(IR)+l RDMAT 60 RI(NZ)=IR RDMAT 61 IF(IC-I-1) 150,140,120 RDMAT 62 120 IS=I+1 RDMAT 63 IE=IC-1 RDMAT 64 DO 130 J=IS,IE RDMAT 65 130 CS(J)=NZ RDMAT 66 140 CS(IC)=NZ RDMAT 67 I=IC RDMAT 68 150 CONTINUE RDMAT 69 PRINT 910 RDMAT 70 910 FORMAT(1X,9(1H#)" TOO MANY NONZEROS...FIX DIMENSIONS") RDMAT 71 STOP RDMAT 72 C RDMAT 73 800 PRINT 920,M RDMAT 74 920 FORMAT(1X,9(1H)" THE MATRIX SIZE IS"I6"...FIX DIMENSIONS") RDMAT 75 STOP RDMAT 76 C RDMAT 77 180 MS=I+1 RDMAT 78 DO 190 J=MS,M1 RDMAT 79 190 CS(J)=NZ RDMAT 80 NZ=NZ-1 RDMAT 81 C RDMAT 82 RETURN RDMAT 83 END RDMAT 84 B.3 SUBROUTINE COLROW (M,CS,RI,RC,RS,CI) COLROW 2 C COLROW 3 C PURPOSE COLROW 4 C GIVEN THE COLUMNWISE STORAGE OF A SPARSE MATRIX, COLR-)W 5 C COLROW CREATES ROW WISE REPRESENTATION ARRAYS (OR VICE VERSA). COLROW 6 C THE NEW ARRAYS ARE ADDITIONAL AND DO NOT REPLACE THE OLD ONES. COLROW 7 C COLROW C COLROW 9 INTEGER CS(M),RI(1),RC(M),RS(M),CI(1) COLROW 10 C COLROW 11 C INPUT ARGUMENTS COLROW 12 C M SIZE OF THE MATRIX COLROW 13 C CS(M+1) COLUMNS STARTING POINTS IN RI( ). COLROW 14 C RI( ) COLUMN WISE LIST OF ROW INDICES. COLROW 15 C RC(M) NONZEROS COUNT FOR ROWS. COLROW 16 C COLROW 17 C OUTPUT ARGUMENTS COLROW 18 C RS(M+1) ROWS STARTING POINTS IN CI( ). COLROW 19 C CI( ) ROW WISE LIST OF COLUMN INDICES. COLROW 20 C COLROW 21 C COLROW 22 M1=M+1 COLROW 23 RS(1)=1 COL.ROW 24 RS(2)=1 COLROW 25 DO 10 I=3,M1 "OLROW 26 10 RS(I)=RS(I-1)+RC(I-2) COLROW 27 C COLROW 28 DO 100 J=1,M COLROW 29 IF(CS(J).EQ.CS(J+1)) GOTO 100 COLROW 30 KS=CS(J) COLROW 31 KE=CS(J+1)-1 COLROW 32 DO 20 K=KS,KE COLROW 33 I=RI(K)+1 COLROW 34 KK=RS(I) COLROW 35 CI(KK)=J COLROW 36 20 RS(I)=RS(I)+1 COLROW 37 100 CONTINUE COLROW 38 C COLROW 39 RETURN COLROW 40 END COLROW 41 B.4 SUBROUTINE PRSYM(CI,RS,IRS,IRE,ICS,ICE,MX,RNTO,CNTO,COTN) PRSYM 2 C PRSYM 3 C PRSYM PRINTS A PICTURE OF A MATRIX, OR ANY SUBMATRIX OF IT, PRSYM 4 C WITH THE NONZERO ELEMENTS MARKED. EACH SET OF 125 COLUMNS 15 PRSYM 5 C PRINTED SEPERATELY. ROWS AND COLUMNZ MAY BE PERMUTED. PRSYM 6 C PRSYM 7 INTEGER CI(1),RS(1),RNTO(1),CNTO(1),COTN(1) PRSYM 8 C PRSYM 9 C INPUT ARGUMENTS (REQUIRED DIMENSION) PRSYM 110 C CI( ) COLUMN INDICES OF NONZERO ELEMENTS, GIVEN nOW WISE. PRSYM 11 C RS(NR+1) PCINTERS TO THE BEGINNING OF ROWS IN CI( ). PRSYM 12 C IRS/ICS FIRST ROW/COLUMN POSITION IN THE PICTURE PRSYM 13 C IRE/ICE LAST ROW/COLUMN POSITION PRSYM 14 C MX MAXIMUM COLUMN INDEX TO BE PRINTED PRSYM 15 C RNTO( ) ROWS INVERSE PERMUTATION PRSYM 16 C CNTO( ) COLUMNS INVERSE PERMUTATION PRSYM 17 C COTN( ) COLUMNS PERMUTATION PRSYM 18 C PRSYM 19 c PRSYM 20 DIMENSION IW(125),ICHAR(11) PRSYM 21 DATA ICHAR/1HO,1H1,1H2,1H3,1H4,1H5,1H5,1H7,1H8,1H9,1H / PRSYM 22 C PRSYM 23 C FIND THE NUMBER OF DECIMAL PLACES REQUIRED TO REPRESENT MX PRSYM 24 C PRSYM 25 IQUOT=MX PRSYM 26 DO 110 I=1,5 PRSYM 27 IQUOT=IQUOT/10 PRSYM 28 IF(IQUCT .LE. 9) GO TO 120 PRSYM 29 110 CONTINUE PRSYM 30 120 NDGT=I+1 PRSYM 31 IF(IQUOT .LE. 0) NDGT=1 PRSYM 32 C PRSYM 33 C PRINT COLUMN INDICES PRSYM 34 C PRSYM 35 JS=ICS PRSYM 36 JE=JS+124 PRSYM 37 130 IF (ICE.LT.JE) JE=ICE PRSYM 3 WRITE (6,910) JS,JE PRSYM 39 910 FORMAT(1H1,1X,7HCOLUMNS,I4,3H TO,14 // ) PRSYM 40 NN=100*(NDGT+1) PRSYM 41 DO 150 K=1,NDGT PRSYM 42 NN=NN/10 PRSYM 43 NNN=NN/10 ?RSYM 44 L=0 ?RSYM 45 DO 14C J=JS,JE PRSYM 46 L=L+1 PRSYM 47 ID=MOD(CNTO(J),,NN)/NNN PRSYM 46 IW(L)=ICHAR(ID+1) PRSYM 49 IF(CNTO(J)/NN .EQ. 0 .AND. ID .EQ. 0) IW(L)=lCHAR(11) PRSYM 50 140 CONTINUE PRSYM 51 WRITE (6,920) (IW(I),I=1,L) ?RSYM 52 920 FORMAT(5X,125Al) PRSYM 53 150 CONTINUE PRSYM 54 C PRSYM 55 C PRINT THE MATRIX RCW BY ROW PRSYM 56 C PRSYM 57 M=5 PRSYM 58 B.5 DO 200 I=IRS,IRE PRSYM 59 M=M+1 PRSYM 60 DO 160 K=1,L PRSYM 61 160 IW(K)=1H PRSYM 62 JMP=5 PRSYM 63 IF (M.LE.5.AND.I.NE.IRE) GOTO 170 PRSYM 64 M=1 PRSYM 65 JMP=1 PRSYM 66 170 DO 180 K=1,L,JMP PRSYM 67 180 IW(K)=1H. PRSYM 68 IF (JE.EQ.ICE) IW(L)=1H. PRSYM 69 C PRSYM 70 IR=RNTO(I) PRSYM 71 KS=RS(IR) PRSYM 72 KE=RS(IR+1)-1 PRSYM 73 C PRSYM 74 DO 190 K=KS,KE PRSYM 75 LL=CI(K) PRSYM 76 KK=COTN(LL) PRSYM 77 IF(KK .LT. JS .OR. KK .GT. JE) GO TO 190 PRSYM 78 KK=KK-JS+1 PRSYM 79 IW(KK)=1HX PRSYM 80 190 CONTINUE PRSYM 81 WRITE (6,930) RNTO(I),(IW(K),K=1,L) PRSYM 82 930 FORMAT(I4,1X,125Al) PRSYM 83 200 CONTINUE PRSYM 84 C PRSYM 85 IF (JE.EQ.ICE) RETURN PRSYM 86 JS=JE+1 PRSYM 87 JE=JE+125 PRSYM 88 GOTO 130 PRSYM 89 C PRSYM 90 END PRSYM 91