1 : C GPSKC SOURCE CHAT 05/01/13 00:19:14 5004 2 : C ALGORITHM 582, COLLECTED ALGORITHMS FROM ACM. 3 : C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.8, NO. 2, 4 : C JUN., 1982, P. 190. 5 : C ============================================================== 6 : C 7 : C GIBBS-POOLE-STOCKMEYER AND GIBBS-KING ALGORITHMS ... 8 : C 9 : C 1. SUBROUTINES GPSKCA, GPSKCB, ..., GPSKCQ WHICH IMPLEMENT 10 : C GIBBS-POOLE-STOCKMEYER AND GIBBS-KING ALGORITHMS ... 11 : C 12 : C 2. SAMPLE DRIVER PROGRAM 13 : C 14 : C 3. SAMPLE TEST PROBLEMS 15 : C 16 : C 4. OUTPUT PRODUCED BY SAMPLE DRIVER ON SAMPLE TEST PROBLEMS 17 : C 18 : C ALL OF THE ABOVE ARE IN 80 COLUMN FORMAT. THE FIRST TWO 19 : C SECTIONS HAVE SEQUENCE NUMBERS IN COLUMNS 73 TO 80. THE 20 : C THIRD SECTION HAS DATA IN ALL 80 COLUMNS. THE LAST SECTION 21 : C HAS CARRIAGE CONTROL CHARACTERS IN COLUMN 1 AND DATA IN ALL 22 : C 80 COLUMNS. 23 : C 24 : C THESE FOUR SECTIONS OF THE FILE ARE SEPARATED BY SINGLE CARDS 25 : C OF THE FORM 'C === SEPARATOR ===' IN COLUMNS 1 TO 19 26 : C 27 : C ============================================================== 28 : C 29 : C THE SAMPLE TEST PROBLEMS INCLUDED WITH THIS CODE PROVIDE A 30 : C MINIMAL CHECKOUT OF THE FUNCTIONING OF THE CODE. THE TEST 31 : C PROBLEMS HAVE NUMERICAL VALUES WITH THE INTEGER PART BEING 32 : C THE ROW INDEX, THE FRACTIONAL PART THE COLUMN INDEX, OF THE 33 : C MATRIX AFTER GPS(K) REORDERING. THE TEST OUTPUT INCLUDES A 34 : C LISTING OF THE REORDERED MATRIX, IN WHICH THE NUMERIC VALUES 35 : C SHOULD APPEAR IN CORRECT POSITIONS, INTERMINGLED WITH ZEROES 36 : C WHICH REPRESENT FILL IN THE SPARSE MATRIX FACTORIZATION. 37 : C 38 : C ============================================================== 39 : C 40 : C === SEPARATOR === BEGINNING OF GPS AND GK ALGORITHMS 41 : SUBROUTINE GPSKCA (N, DEGREE, RSTART, CONNEC, OPTPRO, WRKLEN, 42 : 1 PERMUT, WORK, BANDWD, PROFIL, ERROR, SPACE) 43 : C 44 : C ================================================================== 45 : C ================================================================== 46 : C = = 47 : C = B A N D W I D T H OR P R O F I L E R E D U C T I O N = 48 : C = FOR A SPARSE AND (STRUCTURALLY) SYMMETRIC MATRIX, = 49 : C = USING EITHER = 50 : C = = 51 : C = THE GIBBS-POOLE-STOCKMEYER ALGORITHM (BANDWIDTH REDUCTION) = 52 : C = OR = 53 : C = THE GIBBS-KING ALGORITHM (PROFILE REDUCTION) = 54 : C = = 55 : C ================================================================== 56 : C ================================================================== 57 : C = THIS CODE SUPERSEDES TOMS ALGORITHMS 508 AND 509 IN THE = 58 : C = COLLECTED ALGORITHMS OF THE ACM (CALGO). = 59 : C ================================================================== 60 : C ================================================================== 61 : C 62 : C ------------------- 63 : C P A R A M E T E R S 64 : C ------------------- 65 : C 66 : IMPLICIT INTEGER(I-N) 67 : INTEGER N, RSTART(N), WRKLEN, BANDWD, PROFIL, ERROR, SPACE 68 : C 69 : CIBM INTEGER *2 DEGREE(N), CONNEC(1), PERMUT(N), WORK(WRKLEN) 70 : INTEGER DEGREE(N), CONNEC(1), PERMUT(N), WORK(WRKLEN) 71 : C 72 : LOGICAL OPTPRO 73 : C 74 : C ------------------------------------------------------------------ 75 : C 76 : C INPUT PARAMETERS: 77 : C ----- ---------- 78 : C 79 : C N -- THE DIMENSION OF THE MATRIX 80 : C 81 : C DEGREE, 82 : C RSTART, 83 : C CONNEC -- DESCRIBE THE STRUCTURE OF THE SPARSE MATRIX. 84 : C DEGREE(I) SPECIFIES THE NUMBER OF NON-ZERO 85 : C OFF-DIAGONAL ENTRIES IN THE I-TH ROW OF THE 86 : C SPARSE MATRIX. THE COLUMN INDICES OF THESE 87 : C ENTRIES ARE GIVEN IN CONSECUTIVE LOCATIONS IN 88 : C CONNEC, STARTING AT LOCATION RSTART(I). 89 : C IN OTHER WORDS, THE INDICES OF THE NON-ZERO 90 : C OFF-DIAGONAL ELEMENTS OF THE I-TH ROW ARE FOUND 91 : C IN: 92 : C CONNEC (RSTART(I)), 93 : C CONNEC (RSTART(I) + 1), 94 : C . . . 95 : C CONNEC (RSTART(I) + DEGREE(I) - 1) 96 : C 97 : C DIMENSIONS: 98 : C RSTART IS DIMENSION N (OR LONGER). 99 : C DEGREE IS DIMENSION N (OR LONGER). 100 : C CONNEC IS DIMENSION ROUGHLY THE NUMBER OF NON- 101 : C ZERO ENTRIES IN THE MATRIX. 102 : C 103 : C OPTPRO -- .TRUE. IF REDUCING THE PROFILE OF THE MATRIX 104 : C IS MORE IMPORTANT THAN REDUCING THE 105 : C BANDWIDTH 106 : C .FALSE. IF BANDWIDTH REDUCTION IS MOST IMPORTANT 107 : C 108 : C WRKLEN -- THE ACTUAL LENGTH OF THE VECTOR WORK AS SUPPLIED 109 : C BY THE USER. SEE THE DISCUSSION OF THE WORKSPACE 110 : C 'WORK' BELOW FOR TYPICAL STORAGE REQUIREMENTS. 111 : C THE VALUE OF WRKLEN WILL BE USED TO ENSURE THAT 112 : C THE ROUTINE WILL NOT USE MORE STORAGE THAN IS 113 : C AVAILABLE. IF NOT ENOUGH SPACE IS GIVEN IN WORK 114 : C TO PERMIT A SOLUTION TO BE FOUND, THE ERROR FLAG 115 : C WILL BE SET AND FURTHER COMPUTATION STOPPED. 116 : C 117 : C 118 : C INPUT AND OUTPUT PARAMETER: 119 : C ----- --- ------ --------- 120 : C 121 : C PERMUT -- ON INPUT, AN ALTERNATIVE REORDERING FOR THE 122 : C ROWS AND COLUMNS OF THE MATRIX. PERMUT(I) GIVES 123 : C THE POSITION IN WHICH ROW AND COLUMN I SHOULD 124 : C BE PLACED TO REDUCE THE BANDWIDTH OR THE PROFILE. 125 : C IF THE USER HAS NO ALTERNATIVE TO THE NATURAL 126 : C HE SHOULD INITIALIZE PERMUT TO BE THE IDENTITY 127 : C PERMUTATION PERMUT(I) = I . 128 : C 129 : C ON OUTPUT, PERMUT WILL CONTAIN THE PERMUTATION 130 : C FOR REORDERING THE ROWS AND COLUMNS WHICH REDUCES 131 : C THE BANDWIDTH AND/OR PROFILE. THE RESULT WILL BE 132 : C THE REORDERING FOUND BY 'GPSKCA' OR THE REORDERING 133 : C GIVEN BY THE USER IN 'PERMUT', WHICHEVER DOES THE 134 : C JOB BETTER. 135 : C 136 : C 137 : C OUTPUT PARAMETERS: 138 : C ------ ---------- 139 : C 140 : C WORK -- A TEMPORARY STORAGE VECTOR, OF LENGTH SOMEWHAT 141 : C GREATER THAN 3N. THE SPACE BEYOND 3N REQUIRED 142 : C IS PROBLEM-DEPENDENT. ANY PROBLEM CAN BE SOLVED 143 : C IN 6N+3 LOCATIONS. 144 : C MOST PROBLEMS CAN BE REORDERED WITH 4N LOCATIONS 145 : C IN 'WORK'. IF SPACE IS NOT A CONSTRAINT, PROVIDE 146 : C 6N+3 LOCATIONS IN 'WORK'. OTHERWISE, PROVIDE AS 147 : C MUCH MORE THAN 3N AS IS CONVENIENT AND CHECK THE 148 : C ERROR FLAG AND SPACE REQUIRED PARAMETERS (SEE BELOW) 149 : C 150 : C ON OUTPUT, THE 1ST N LOCATIONS OF WORK WILL BE 151 : C A LISTING OF THE ORIGINAL ROW AND COLUMN INDICES AS 152 : C THEY APPEAR IN THE COMPUTED REORDERING. 153 : C LOCATIONS N+1, ... , 2N OF WORK WILL CONTAIN 154 : C THE NEW POSITIONS FOR THE EQUATIONS IN THE ORDER 155 : C FOUND BY GPSKCA. THUS, THE TWO VECTORS ARE INVERSE 156 : C PERMUTATIONS OF EACH OTHER. IF THE ORDERING 157 : C FOUND BY THIS ALGORITHM IS BETTER THAN THE USER- 158 : C SUPPLIED ORDER, THE SECOND PERMUTATION VECTOR IS 159 : C IDENTICAL TO THE RESULT RETURNED IN 'PERMUT'. 160 : C 161 : C BANDWD -- THE BANDWIDTH OF THE MATRIX WHEN ROWS AND COLUMNS 162 : C ARE REORDERED IN THE ORDERING RETURNED IN PERMUT. 163 : C 164 : C PROFIL -- THE PROFILE OF THE MATRIX WHEN ROWS AND COLUMNS ARE 165 : C REORDERED IN THE ORDERING RETURNED IN PERMUT. 166 : C 167 : C ERROR -- WILL BE EQUAL TO ZERO IF A NEW NUMBERING COULD BE 168 : C FOUND IN THE SPACE PROVIDED. OTHERWISE, ERROR 169 : C WILL BE SET TO A POSITIVE ERROR CODE (SEE TABLE 170 : C GIVEN BELOW). IF THE REORDERING ALGORITHM HAS BEEN 171 : C STOPPED BY LACK OF WORKSPACE, THE SPACE PARAMETER 172 : C WILL BE SET TO THE NUMBER OF ADDITIONAL LOCATIONS 173 : C REQUIRED TO COMPLETE AT LEAST THE NEXT PHASE OF 174 : C THE ALGORITHM. 175 : C 176 : C WHENEVER A NON-ZERO VALUE FOR ERROR IS GIVEN 177 : C PERMUT WILL RETAIN THE VALUES PROVIDED BY THE USER 178 : C AND THE SCALARS BANDWD AND PROFIL WILL BE SET TO 179 : C OUTRAGEOUS VALUES. IT IS THE USER'S RESPONSIBILITY 180 : C TO CHECK THE STATUS OF ERROR. 181 : C 182 : C SPACE -- WILL INDICATE EITHER HOW MUCH SPACE THE REORDERING 183 : C ACTUALLY REQUIRED OR HOW MUCH SPACE WILL BE 184 : C REQUIRED TO COMPLETE THE NEXT PHASE OF THE 185 : C REORDERING ALGORITHM. THE POSSIBLE OUTCOMES ARE .. 186 : C 187 : C ERROR = 0 SPACE IS THE MINIMAL VALUE FOR 188 : C WRKLEN REQUIRED TO REORDER 189 : C THIS MATRIX AGAIN. 190 : C 191 : C ERROR <> 0 SPACE IS THE MINIMUM NUMBER 192 : C DUE TO LACK OF OF EXTRA WORKSPACE REQUIRED 193 : C WORKSPACE TO CONTINUE THE REORDERING 194 : C ALGORITHM ON THIS MATRIX. 195 : C 196 : C ERROR <> 0 SPACE = -1 197 : C DUE TO ERROR 198 : C IN DATA STRUCTURES 199 : C 200 : C 201 : C ================================================================== 202 : C 203 : C ---------------------- 204 : C E R R O R C O D E S 205 : C ---------------------- 206 : C 207 : C ERROR CODES HAVE THE FORM 0XY OR 1XY. 208 : C 209 : C ERRORS OF THE FORM 1XY RESULT FROM INADEQUATE WORKSPACE. 210 : C 211 : C ERRORS OF THE FORM 0XY ARE INTERNAL PROGRAM CHECKS, WHICH 212 : C MOST LIKELY OCCUR BECAUSE THE CONNECTIVITY STRUCTURE OF THE 213 : C MATRIX IS REPRESENTED INCORRECTLY (E.G., THE DEGREE OF 214 : C A NODE IS NOT CORRECT OR NODE I IS CONNECTED TO NODE J, 215 : C BUT NOT CONVERSELY). 216 : C 217 : C THE LAST DIGIT (Y) IS MAINLY USEFUL FOR DEBUGGING THE 218 : C THE REORDERING ALGORITHM. THE MIDDLE DIGIT (X) INDICATES 219 : C HOW MUCH OF THE ALGORITHM HAS BEEN PERFORMED. 220 : C THE TABLE BELOW GIVES THE CORRESPONDENCE BETWEEN THE 221 : C VALUES OF X AND THE STRUCTURE OF THE ALGORITHM. 222 : C X = 0 INITIAL PROCESSING 223 : C X = 1 COMPUTING PSEUDO-DIAMETER (ALGORITHM I) 224 : C X = 2 TRANSITION BETWEEN ALGORITHM I AND II 225 : C X = 3 COMBINING LEVEL STRUCTURES (ALGORITHM II) 226 : C X = 4 TRANSITION BETWEEN ALGORITHM II AND III 227 : C X = 5 BANDWIDTH NUMBERING (ALGORITHM IIIA) 228 : C X = 6 PROFILE NUMBERING (ALGORITHM IIIB) 229 : C X = 7 FINAL BANDWIDTH/PROFILE COMPUTATION 230 : C 231 : C ================================================================== 232 : C 233 : C --------------------- --------------- 234 : C A L T E R N A T I V E V E R S I O N S 235 : C --------------------- --------------- 236 : C 237 : C SHORT INTEGER VERSION 238 : C 239 : C ON MACHINES WITH TWO OR MORE PRECISIONS FOR INTEGERS, 240 : C ALL OF THE INPUT ARRAYS EXCEPT 'RSTART' CAN BE CONVERTED 241 : C TO THE SHORTER PRECISION AS LONG AS THAT SHORTER PRECISION 242 : C ALLOWS NUMBERS AS LARGE AS 'N'. A VERSION OF THIS CODE 243 : C SUITABLE FOR USE ON IBM COMPUTERS (INTEGER * 2) IS EMBEDDED 244 : C AS COMMENTS IN THIS CODE. ALL SUCH COMMENTS HAVE THE 245 : C CHARACTERS 'CIBM' IN THE FIRST FOUR COLUMNS, AND PRECEDE THE 246 : C EQUIVALENT STANDARD CODE WHICH THEY WOULD REPLACE. 247 : C 248 : C CONNECTIVITY COMPATIBILITY VERSION 249 : C 250 : C THE ORIGINAL (1976) TOMS CODE 'REDUCE' USED A LESS STORAGE 251 : C EFFICIENT FORMAT FOR THE CONNECTIVITY TABLE 'CONNEC'. 252 : C THE 1976 CODE USED A RECTANGULAR MATRIX OF DIMENSIONS 253 : C N BY MAXDGR, WHERE MAXDGR IS AT LEAST AS LARGE AS 254 : C THE MAXIMUM DEGREE OF ANY NODE IN THE GRAPH OF THE MATRIX. 255 : C THE FORMAT USED IN THE CURRENT CODE IS OFTEN SUBSTANTIALLY 256 : C MORE EFFICIENT. HOWEVER, FOR USERS FOR WHOM CONVERSION WILL 257 : C BE DIFFICULT OR IMPOSSIBLE, TWO ALTERNATIVES ARE .. 258 : C 1. SIMPLY NOTE THAT CHANGING THE ORDER OF SUBSCRIPTS 259 : C IN A RECTANGULAR CONNECTION TABLE WILL ENABLE YOU 260 : C TO USE THE NEW VERSION. THIS SUBROUTINE WILL ACCEPT A 261 : C RECTANGULAR CONNECTION TABLE OF DIMENSIONS 262 : C MAXDGR BY N, 263 : C PROVIDED THAT RSTART(I) IS SET TO (I-1)*MAXDGR + 1. 264 : C 2. THE AUTHOR WILL MAKE AVAILABLE A VARIANT VERSION 265 : C 'GPSKRA', WHICH EXPECTS THE ADJACENCY MATRIX OR 266 : C CONNECTIVITY TABLE IN THE SAME FORM AS DID 'REDUCE'. 267 : C THIS VERSION CAN BE OBTAINED BY WRITING TO .. 268 : C JOHN GREGG LEWIS 269 : C BOEING COMPUTER SERVICES COMPANY 270 : C MAIL STOP 9C-01 271 : C P.O. BOX 24346 272 : C SEATTLE, WA 98124 273 : C PLEASE INCLUDE A DESCRIPTION OF THE COMPUTING 274 : C ENVIRONMENT ON WHICH YOU WILL BE USING THE CODE. 275 : C 276 : C ================================================================== 277 : C 278 : INTEGER I, INC1, INC2, AVAIL, NXTNUM, LOWDG, STNODE, NLEFT, 279 : 1 TREE1, TREE2, DEPTH, EMPTY, STOTAL, REQD, CSPACE, 280 : 2 LVLLST, LVLPTR, ACTIVE, RVNODE, WIDTH1, WIDTH2, MXDG 281 : C 282 : LOGICAL REVRS1, ONEIS1 283 : C 284 : C ================================================================== 285 : C 286 : C << NUMBER ANY DEGREE ZERO NODES >> 287 : C 288 : C WHILE << SOME NODES YET UNNUMBERED >> DO 289 : C << FIND A PSEUDO-DIAMETER OF THE MATRIX GRAPH >> 290 : C << CONVERT FORM OF LEVEL TREES >> 291 : C << COMBINE LEVEL TREES INTO ONE LEVEL STRUCTURE >> 292 : C << CONVERT FORM OF LEVEL STRUCTURE >> 293 : C IF OPTPRO THEN 294 : C << RENUMBER BY KING ALGORITHM >> 295 : C ELSE 296 : C << RENUMBER BY REVERSE CUTHILL-MCKEE ALGORITHM >> 297 : C 298 : C ================================================================== 299 : C 300 : C ... INITIALIZE COUNTERS, THEN NUMBER ANY NODES OF DEGREE 0. 301 : C THE LIST OF NODES, BY NEW NUMBER, WILL BE BUILT IN PLACE AT 302 : C THE FRONT OF THE WORK AREA. 303 : C 304 : NXTNUM = 1 305 : ERROR = 0 306 : SPACE = 2*N 307 : C 308 : MXDG = 0 309 : DO 300 I = 1, N 310 : IF (DEGREE(I)) 6000, 100, 200 311 : 100 WORK(NXTNUM) = I 312 : NXTNUM = NXTNUM + 1 313 : GO TO 300 314 : 200 IF (DEGREE(I) .GT. MXDG) MXDG = DEGREE(I) 315 : 300 CONTINUE 316 : C 317 : C 318 : C ============================== 319 : C ... WHILE NXTNUM <= N DO ... 320 : C ============================== 321 : C 322 : 1000 IF ( NXTNUM .GT. N ) GO TO 2000 323 : C 324 : C ... FIND AN UNNUMBERED NODE OF MINIMAL DEGREE 325 : C 326 : LOWDG = MXDG + 1 327 : STNODE = 0 328 : DO 400 I = 1, N 329 : IF ( (DEGREE(I) .LE. 0) .OR. (DEGREE(I) .GE. LOWDG) ) 330 : 1 GO TO 400 331 : LOWDG = DEGREE(I) 332 : STNODE = I 333 : 400 CONTINUE 334 : C 335 : IF ( STNODE .EQ. 0 ) GO TO 6100 336 : C 337 : C ... SET UP POINTERS FOR THREE LISTS IN WORK AREA, THEN LOOK 338 : C FOR PSEUDO-DIAMETER, BEGINNING WITH STNODE. 339 : C 340 : AVAIL = (WRKLEN - NXTNUM + 1) / 3 341 : NLEFT = N - NXTNUM + 1 342 : SPACE = MAX0 (SPACE, NXTNUM + 3*N - 1) 343 : IF ( AVAIL .LT. N ) GO TO 5200 344 : C 344 : CALL GPSKCB (N, DEGREE, RSTART, CONNEC, AVAIL, NLEFT, 345 : 1 STNODE, RVNODE, WORK(NXTNUM), TREE1, TREE2, 346 : 2 ACTIVE, DEPTH, WIDTH1, WIDTH2, 347 : 3 ERROR, SPACE) 348 : IF ( ERROR .NE. 0 ) GO TO 5000 349 : SPACE = MAX0 (SPACE, NXTNUM + 3*(ACTIVE+DEPTH+1) - 1) 350 : C 351 : C ... DYNAMIC SPACE CHECK FOR MOST OF REMAINDER OF ALGORITHM 352 : C 353 : REQD = MAX0 (NXTNUM + 2*N + 3*DEPTH - 1, 3*N + 2*DEPTH + 1) 354 : SPACE = MAX0 (SPACE, REQD) 355 : IF ( WRKLEN .LT. REQD ) GO TO 5300 356 : C 357 : C 358 : C ... OUTPUT FROM GPSKCB IS A PAIR OF LEVEL TREES, IN THE FORM 359 : C OF LISTS OF NODES BY LEVEL. CONVERT THIS TO TWO LISTS OF 360 : C OF LEVEL NUMBER BY NODE. AT THE SAME TIME PACK 361 : C STORAGE SO THAT ONE OF THE LEVEL TREE VECTORS IS AT THE 362 : C BACK END OF THE WORK AREA. 363 : C 364 : LVLPTR = NXTNUM + AVAIL - DEPTH 364 : CALL GPSKCE (N, AVAIL, ACTIVE, DEPTH, WRKLEN, WORK(NXTNUM), 365 : 1 WORK(LVLPTR), WORK(1), NXTNUM, TREE1, 366 : 2 TREE2, WIDTH1, WIDTH2, ONEIS1, ERROR, SPACE) 367 : IF ( ERROR .NE. 0 ) GO TO 5000 368 : IF (( TREE1 .NE. WRKLEN - N + 1 ) .OR. (TREE2 .NE. NXTNUM)) 369 : 1 GO TO 6200 370 : C 371 : C ... COMBINE THE TWO LEVEL TREES INTO A MORE GENERAL 372 : C LEVEL STRUCTURE. 373 : C 374 : AVAIL = WRKLEN - NXTNUM + 1 - 2*N - 3*DEPTH 375 : STOTAL = N + NXTNUM 376 : EMPTY = STOTAL + DEPTH 377 : INC1 = TREE1 - DEPTH 378 : INC2 = INC1 - DEPTH 379 : C 379 : CALL GPSKCG (N, DEGREE, RSTART, CONNEC, ACTIVE, WIDTH1, 380 : 1 WIDTH2, WORK(TREE1), WORK(TREE2), WORK(EMPTY), 381 : 2 AVAIL, DEPTH, WORK(INC1), WORK(INC2), 382 : 3 WORK(STOTAL), ONEIS1, REVRS1, ERROR, CSPACE) 383 : C 384 : IF ( ERROR .NE. 0 ) GO TO 5000 385 : SPACE = MAX0 (SPACE, NXTNUM + CSPACE - 1) 386 : C 387 : C ... COMBINED LEVEL STRUCTURE IS REPRESENTED BY GPSKCG AS 388 : C A VECTOR OF LEVEL NUMBERS. FOR RENUMBERING PHASE, 389 : C CONVERT THIS ALSO TO THE INVERSE PERMUTATION. 390 : C 391 : LVLPTR = TREE1 - (DEPTH + 1) 392 : LVLLST = LVLPTR - ACTIVE 393 : IF ( STOTAL + DEPTH .GT. LVLPTR ) GO TO 6300 394 : C 394 : CALL GPSKCI (N, ACTIVE, DEPTH, WORK(TREE1), WORK(LVLLST), 395 : 1 WORK(LVLPTR), WORK(STOTAL), ERROR, SPACE) 396 : IF (ERROR .NE. 0) GO TO 5000 397 : C 398 : C ... NOW RENUMBER ALL MEMBERS OF THIS COMPONENT USING 399 : C EITHER A REVERSE CUTHILL-MCKEE OR A KING STRATEGY, 400 : C AS PROFILE OR BANDWIDTH REDUCTION IS MORE IMPORTANT. 401 : C 402 : IF ( OPTPRO ) GO TO 500 402 : CALL GPSKCJ (N, DEGREE, RSTART, CONNEC, ACTIVE, 403 : 1 WORK(NXTNUM), STNODE, RVNODE, REVRS1, DEPTH, 404 : 2 WORK(LVLLST), WORK(LVLPTR), WORK(TREE1), 405 : 3 ERROR, SPACE) 406 : IF ( ERROR .NE. 0 ) GO TO 5000 407 : NXTNUM = NXTNUM + ACTIVE 408 : GO TO 600 409 : C 409 : 500 CALL GPSKCK (N, DEGREE, RSTART, CONNEC, LVLLST-1, NXTNUM, 410 : 1 WORK, ACTIVE, DEPTH, WORK(LVLLST), 411 : 2 WORK(LVLPTR), WORK(TREE1), ERROR, SPACE) 412 : IF ( ERROR .NE. 0 ) GO TO 5000 413 : C 414 : C ========================================================= 415 : C ... END OF WHILE LOOP ... REPEAT IF GRAPH IS DISCONNECTED 416 : C ========================================================= 417 : C 418 : 600 GO TO 1000 419 : C 420 : C ... CHECK WHETHER INITIAL NUMBERING OR FINAL NUMBERING 421 : C PROVIDES BETTER RESULTS 422 : C 423 : 2000 IF (WRKLEN .LT. 2*N) GO TO 5400 424 : C 425 : IF (OPTPRO) GO TO 2100 425 : CALL GPSKCL (N, DEGREE, RSTART, CONNEC, WORK(1), WORK(N+1), 426 : 1 PERMUT, BANDWD, PROFIL, ERROR, SPACE) 427 : GO TO 2200 428 : C 428 : 2100 CALL GPSKCM (N, DEGREE, RSTART, CONNEC, WORK(1), WORK(N+1), 429 : 1 PERMUT, BANDWD, PROFIL, ERROR, SPACE) 430 : C 431 : 2200 RETURN 432 : C 433 : C 434 : C . . . E R R O R D I A G N O S T I C S 435 : C --------------------------------- 436 : C 437 : C ... ERROR DETECTED BY LOWER LEVEL ROUTINE. MAKE SURE THAT SIGNS 438 : C OF DEGREE ARE PROPERLY SET 439 : C 440 : 5000 DO 5100 I = 1, N 441 : IF (DEGREE(I) .LT. 0) DEGREE(I) = -DEGREE(I) 442 : 5100 CONTINUE 443 : C 444 : BANDWD = -1 445 : PROFIL = -1 446 : RETURN 447 : C 448 : C ... STORAGE ALLOCATION ERRORS DETECTED IN THIS ROUTINE 449 : C 450 : 5200 ERROR = 101 451 : SPACE = -1 452 : GO TO 5000 453 : C 454 : 5300 ERROR = 102 455 : SPACE = -1 456 : GO TO 5000 457 : C 458 : 5400 ERROR = 10 459 : SPACE = 2*N - WRKLEN 460 : GO TO 5000 461 : C 462 : C ... DATA STRUCTURE ERRORS DETECTED IN THIS ROUTINE 463 : C 464 : 6000 ERROR = 1 465 : GO TO 6900 466 : C 467 : 6100 ERROR = 2 468 : GO TO 6900 469 : C 470 : 6200 ERROR = 3 471 : GO TO 6900 472 : C 473 : 6300 ERROR = 4 474 : C 475 : 6900 SPACE = -1 476 : GO TO 5000 477 : END 478 : SUBROUTINE GPSKCB (N, DEGREE, RSTART, CONNEC, AVAIL, NLEFT, 479 : 1 STNODE, RVNODE, WORK, FORWD, BESTBK, NNODES, 480 : 2 DEPTH, FWIDTH, BWIDTH, ERROR, SPACE) 481 : C 482 : C ================================================================== 483 : C 484 : C FIND A PSEUDO-DIAMETER OF THE MATRIX GRAPH ... 485 : C 486 : C << BUILD A LEVEL TREE FROM STNODE >> 487 : C REPEAT 488 : C << BUILD A LEVEL TREE FROM EACH NODE 'BKNODE' IN THE 489 : C DEEPEST LEVEL OF STNODE'S TREE >> 490 : C << REPLACE 'STNODE' WITH 'BKNODE' IF A DEEPER AND 491 : C NARROWER TREE WAS FOUND. >> 492 : C UNTIL 493 : C << NO FURTHER IMPROVEMENT MADE >> 494 : C 495 : C ... HEURISTIC ABOVE DIFFERS FROM THE ALGORITHM PUBLISHED IN 496 : C SIAM J. NUMERICAL ANALYSIS, BUT MATCHES THE CODE 497 : C DISTRIBUTED BY TOMS. 498 : C 499 : C 500 : C PARAMETERS : 501 : C 502 : C N, DEGREE, RSTART & CONNEC DESCRIBE THE MATRIX STRUCTURE 503 : C 504 : C WORK -- WORKING SPACE, OF LENGTH 3*AVAIL, USED TO STORE 505 : C THREE LEVEL TREES. 506 : C 507 : C STNODE IS INITIALLY THE NUMBER OF A NODE TO BE USED TO 508 : C START THE PROCESS, TO BE THE ROOT OF THE FIRST TREE. 509 : C ON OUTPUT, STNODE IS THE END OF THE PSEUDO-DIAMETER WHOSE 510 : C LEVEL TREE IS NARROWEST. 511 : C 512 : C RVNODE WILL BE THE OTHER END OF THE PSEUDO-DIAMETER. 513 : C 514 : C NNODES WILL BE THE NUMBER OF NODES IN THIS CONNECTED 515 : C COMPONNENT OF THE MATRIX GRAPH, I.E., THE LENGTH OF 516 : C THE LEVEL TREES. 517 : C 518 : C DEPTH -- THE DEPTH OF THE LEVEL TREES BEING RETURNED, 519 : C I.E., THE LENGTH OF THE PSEUDO-DIAMETER. 520 : C 521 : C ================================================================== 522 : C 523 : C STRUCTURE OF WORKSPACE ... 524 : C 525 : C --------------------------------------------------------------- 526 : C : NUMBERED : TLIST1 PTR1 : TLIST2 PTR2 : TLIST3 PTR3 : 527 : C --------------------------------------------------------------- 528 : C 529 : C TLISTI IS A LIST OF NODES OF LENGTH 'ACTIVE' 530 : C PTRI IS A LIST OF POINTERS INTO TLISTI, OF LENGTH 'DEPTH+1' 531 : C 532 : C ================================================================== 533 : C 534 : IMPLICIT INTEGER(I-N) 535 : INTEGER N, RSTART(N), AVAIL, NLEFT, 536 : 1 STNODE, RVNODE, FORWD, BESTBK, NNODES, DEPTH, FWIDTH, 537 : 4 BWIDTH, ERROR, SPACE 538 : C 539 : CIBM INTEGER *2 DEGREE(N), CONNEC(1), WORK(AVAIL,3) 540 : INTEGER DEGREE(N), CONNEC(1), WORK(AVAIL,3) 541 : C 542 : C ---------------- 543 : C 544 : INTEGER BACKWD, MXDPTH, WIDTH, FDEPTH, LSTLVL, 545 : 1 NLAST, T, I, BKNODE, LSTLVI 546 : C 547 : LOGICAL IMPROV 548 : C 549 : C 550 : C ... BUILD INITIAL LEVEL TREE FROM 'STNODE'. FIND OUT HOW MANY 551 : C NODES LIE IN THE CURRENT CONNECTED COMPONENT. 552 : C 553 : FORWD = 1 554 : BACKWD = 2 555 : BESTBK = 3 556 : C 556 : CALL GPSKCC (N, DEGREE, RSTART, CONNEC, STNODE, AVAIL, NLEFT, 557 : 1 WORK(1,FORWD), NNODES, DEPTH, WIDTH, ERROR, 558 : 2 SPACE) 559 : IF ( ERROR .NE. 0 ) GO TO 5000 560 : C 561 : MXDPTH = AVAIL - NNODES - 1 562 : C 563 : C ========================================== 564 : C REPEAT UNTIL NO DEEPER TREES ARE FOUND ... 565 : C ========================================== 566 : C 567 : 1000 FWIDTH = WIDTH 568 : FDEPTH = DEPTH 569 : LSTLVL = AVAIL - DEPTH + 1 570 : NLAST = WORK (LSTLVL-1, FORWD) - WORK (LSTLVL, FORWD) 571 : LSTLVL = WORK (LSTLVL, FORWD) 572 : BWIDTH = N+1 573 : C 574 : C ... SORT THE DEEPEST LEVEL OF 'FORWD' TREE INTO INCREASING 575 : C ORDER OF NODE DEGREE. 576 : C 576 : CALL GPSKCQ (NLAST, WORK(LSTLVL,FORWD), N, DEGREE, ERROR) 577 : IF (ERROR .NE. 0) GO TO 6000 578 : C 579 : C ... BUILD LEVEL TREE FROM NODES IN 'LSTLVL' UNTIL A DEEPER 580 : C AND NARROWER TREE IS FOUND OR THE LIST IS EXHAUSTED. 581 : C 582 : IMPROV = .FALSE. 583 : DO 1200 I = 1, NLAST 584 : LSTLVI = LSTLVL + I - 1 585 : BKNODE = WORK (LSTLVI, FORWD) 585 : CALL GPSKCD (N, DEGREE, RSTART, CONNEC, BKNODE, AVAIL, 586 : 1 NNODES, MXDPTH, WORK(1,BACKWD), DEPTH, WIDTH, 587 : 2 BWIDTH, ERROR, SPACE) 588 : IF ( ERROR .NE. 0 ) GO TO 5000 589 : C 590 : IF ( DEPTH .LE. FDEPTH ) GO TO 1100 591 : C 592 : C ... NEW DEEPER TREE ... MAKE IT NEW 'FORWD' TREE 593 : C AND BREAK OUT OF 'DO' LOOP. 594 : C 595 : IMPROV = .TRUE. 596 : T = FORWD 597 : FORWD = BACKWD 598 : BACKWD = T 599 : STNODE = BKNODE 600 : GO TO 1300 601 : C 602 : C ... ELSE CHECK FOR NARROWER TREE. 603 : C 604 : 1100 IF ( WIDTH .GE. BWIDTH ) GO TO 1200 605 : T = BESTBK 606 : BESTBK = BACKWD 607 : BACKWD = T 608 : BWIDTH = WIDTH 609 : RVNODE = BKNODE 610 : 1200 CONTINUE 611 : C 612 : C ... END OF REPEAT LOOP 613 : C ---------------------- 614 : C 615 : 1300 IF ( IMPROV ) GO TO 1000 616 : C 617 : DEPTH = FDEPTH 618 : RETURN 619 : C 620 : C ... IN CASE OF ERROR, SIMPLY RETURN ERROR FLAG TO USER. 621 : C 622 : 5000 RETURN 623 : C 624 : 6000 ERROR = 11 625 : SPACE = -1 626 : RETURN 627 : C 628 : END 629 : SUBROUTINE GPSKCC (N, DEGREE, RSTART, CONNEC, STNODE, AVAIL, 630 : 1 NLEFT, LIST, ACTIVE, DEPTH, WIDTH, ERROR, 631 : 2 SPACE) 632 : C 633 : C ================================================================== 634 : C BUILD THE LEVEL TREE ROOTED AT 'STNODE' IN THE SPACE PROVIDED IN 635 : C LIST. CHECK FOR OVERRUN OF SPACE ALLOCATION. 636 : C ================================================================== 637 : C 638 : IMPLICIT INTEGER(I-N) 639 : INTEGER N, RSTART(N), STNODE, AVAIL, NLEFT, 640 : 1 ACTIVE, DEPTH, WIDTH, ERROR, SPACE 641 : C 642 : CIBM INTEGER *2 DEGREE(N), CONNEC(1), LIST(AVAIL) 643 : INTEGER DEGREE(N), CONNEC(1), LIST(AVAIL) 644 : C 645 : C ... PARAMETERS: 646 : C 647 : C INPUT ... 648 : C 649 : C N, DEGREE, RSTART, CONNEC -- DESCRIBE THE MATRIX STRUCTURE 650 : C 651 : C STNODE -- THE ROOT OF THE LEVEL TREE. 652 : C 653 : C AVAIL -- THE LENGTH OF THE WORKING SPACE AVAILABLE 654 : C 655 : C NLEFT -- THE NUMBER OF NODES YET TO BE NUMBERED 656 : C 657 : C LIST -- THE WORKING SPACE. 658 : C 659 : C OUTPUT ... 660 : C 661 : C ACTIVE -- THE NUMBER OF NODES IN THE COMPONENT 662 : C 663 : C DEPTH -- THE DEPTH OF THE LEVEL TREE ROOTED AT STNODE. 664 : C 665 : C WIDTH -- THE WIDTH OF THE LEVEL TREE ROOTED AT STNODE. 666 : C 667 : C ERROR -- ZERO UNLESS STORAGE WAS INSUFFICIENT. 668 : C 669 : C ------------------------------------------------------------------ 670 : C 671 : INTEGER LSTART, NLEVEL, FRONT, J, NEWNOD, PTR, CDGREE, 672 : 1 LFRONT, LISTJ 673 : C 674 : C ... BUILD THE LEVEL TREE USING LIST AS A QUEUE AND LEAVING 675 : C THE NODES IN PLACE. THIS GENERATES THE NODES ORDERED BY LEVEL 676 : C PUT POINTERS TO THE BEGINNING OF EACH LEVEL, BUILDING FROM 677 : C THE BACK OF THE WORK AREA. 678 : C 679 : ACTIVE = 1 680 : DEPTH = 0 681 : WIDTH = 0 682 : ERROR = 0 683 : LSTART = 1 684 : FRONT = 1 685 : LIST (ACTIVE) = STNODE 686 : DEGREE (STNODE) = -DEGREE (STNODE) 687 : LIST (AVAIL) = 1 688 : NLEVEL = AVAIL 689 : C 690 : C ... REPEAT UNTIL QUEUE BECOMES EMPTY OR WE RUN OUT OF SPACE. 691 : C ------------------------------------------------------------ 692 : C 693 : 1000 IF ( FRONT .LT. LSTART ) GO TO 1100 694 : C 695 : C ... FIRST NODE OF LEVEL. UPDATE POINTERS. 696 : C 697 : LSTART = ACTIVE + 1 698 : WIDTH = MAX0 (WIDTH, LSTART - LIST(NLEVEL)) 699 : NLEVEL = NLEVEL - 1 700 : DEPTH = DEPTH + 1 701 : IF ( NLEVEL .LE. ACTIVE ) GO TO 5000 702 : LIST (NLEVEL) = LSTART 703 : C 704 : C ... FIND ALL NEIGHBORS OF CURRENT NODE, ADD THEM TO QUEUE. 705 : C 706 : 1100 LFRONT = LIST (FRONT) 707 : PTR = RSTART (LFRONT) 708 : CDGREE = -DEGREE (LFRONT) 709 : IF (CDGREE .LE. 0) GO TO 6000 710 : DO 1200 J = 1, CDGREE 711 : NEWNOD = CONNEC (PTR) 712 : PTR = PTR + 1 713 : C 714 : C ... ADD TO QUEUE ONLY NODES NOT ALREADY IN QUEUE 715 : C 716 : IF ( DEGREE(NEWNOD) .LE. 0 ) GO TO 1200 717 : DEGREE (NEWNOD) = -DEGREE (NEWNOD) 718 : ACTIVE = ACTIVE + 1 719 : IF ( NLEVEL .LE. ACTIVE ) GO TO 5000 720 : IF ( ACTIVE .GT. NLEFT ) GO TO 6000 721 : LIST (ACTIVE) = NEWNOD 722 : 1200 CONTINUE 723 : FRONT = FRONT + 1 724 : C 725 : C ... IS QUEUE EMPTY? 726 : C ------------------- 727 : C 728 : IF ( FRONT .LE. ACTIVE ) GO TO 1000 729 : C 730 : C ... YES, THE TREE IS BUILT. UNDO OUR MARKINGS. 731 : C 732 : DO 1300 J = 1, ACTIVE 733 : LISTJ = LIST(J) 734 : DEGREE (LISTJ) = -DEGREE (LISTJ) 735 : 1300 CONTINUE 736 : C 737 : RETURN 738 : C 739 : C ... INSUFFICIENT STORAGE ... 740 : C 741 : 5000 SPACE = 3 * ( (NLEFT+1-ACTIVE)*DEPTH / NLEFT + (NLEFT+1-ACTIVE) ) 742 : ERROR = 110 743 : RETURN 744 : C 745 : 6000 ERROR = 12 746 : SPACE = -1 747 : RETURN 748 : C 749 : END 750 : SUBROUTINE GPSKCD (N, DEGREE, RSTART, CONNEC, STNODE, AVAIL, 751 : 1 ACTIVE, MXDPTH, LIST, DEPTH, WIDTH, MAXWID, 752 : 2 ERROR, SPACE) 753 : C 754 : C ================================================================== 755 : C BUILD THE LEVEL TREE ROOTED AT 'STNODE' IN THE SPACE PROVIDED IN 756 : C LIST. OVERFLOW CHECK NEEDED ONLY ON DEPTH OF TREE. 757 : C 758 : C BUILD THE LEVEL TREE TO COMPLETION ONLY IF THE WIDTH OF ALL 759 : C LEVELS IS SMALLER THAN 'MAXWID'. IF A WIDER LEVEL IS FOUND 760 : C TERMINATE THE CONSTRUCTION. 761 : C ================================================================== 762 : C 763 : IMPLICIT INTEGER(I-N) 764 : INTEGER N, RSTART(N), STNODE, AVAIL, ACTIVE, MXDPTH, 765 : 1 DEPTH, WIDTH, MAXWID, ERROR, SPACE 766 : C 767 : CIBM INTEGER *2 DEGREE(N), CONNEC(1), LIST(AVAIL) 768 : INTEGER DEGREE(N), CONNEC(1), LIST(AVAIL) 769 : C 770 : C ... PARAMETERS: 771 : C 772 : C INPUT ... 773 : C 774 : C N, DEGREE, RSTART, CONNEC -- DESCRIBE THE MATRIX STRUCTURE 775 : C 776 : C STNODE -- THE ROOT OF THE LEVEL TREE. 777 : C 778 : C AVAIL -- THE LENGTH OF THE WORKING SPACE AVAILABLE 779 : C 780 : C NLEFT -- THE NUMBER OF NODES YET TO BE NUMBERED 781 : C 782 : C ACTIVE -- THE NUMBER OF NODES IN THE COMPONENT 783 : C 784 : C MXDPTH -- MAXIMUM DEPTH OF LEVEL TREE POSSIBLE IN 785 : C ALLOTTED WORKING SPACE 786 : C 787 : C LIST -- THE WORKING SPACE. 788 : C 789 : C OUTPUT ... 790 : C 791 : C DEPTH -- THE DEPTH OF THE LEVEL TREE ROOTED AT STNODE. 792 : C 793 : C WIDTH -- THE WIDTH OF THE LEVEL TREE ROOTED AT STNODE. 794 : C 795 : C MAXWID -- LIMIT ON WIDTH OF THE TREE. TREE WILL NOT BE 796 : C USED IF WIDTH OF ANY LEVEL IS AS GREAT AS 797 : C MAXWID, SO CONSTRUCTION OF TREE NEED NOT 798 : C CONTINUE IF ANY LEVEL THAT WIDE IS FOUND. 799 : C ERROR -- ZERO UNLESS STORAGE WAS INSUFFICIENT. 800 : C 801 : C ------------------------------------------------------------------ 802 : C 803 : INTEGER LSTART, NLEVEL, FRONT, J, NEWNOD, PTR, BACK, 804 : 1 SPTR, FPTR, LFRONT, LISTJ 805 : C 806 : C ... BUILD THE LEVEL TREE USING LIST AS A QUEUE AND LEAVING 807 : C THE NODES IN PLACE. THIS GENERATES THE NODES ORDERED BY LEVEL 808 : C PUT POINTERS TO THE BEGINNING OF EACH LEVEL, BUILDING FROM 809 : C THE BACK OF THE WORK AREA. 810 : C 811 : BACK = 1 812 : DEPTH = 0 813 : WIDTH = 0 814 : ERROR = 0 815 : LSTART = 1 816 : FRONT = 1 817 : LIST (BACK) = STNODE 818 : DEGREE (STNODE) = -DEGREE (STNODE) 819 : LIST (AVAIL) = 1 820 : NLEVEL = AVAIL 821 : C 822 : C ... REPEAT UNTIL QUEUE BECOMES EMPTY OR WE RUN OUT OF SPACE. 823 : C ------------------------------------------------------------ 824 : C 825 : 1000 IF ( FRONT .LT. LSTART ) GO TO 1100 826 : C 827 : C ... FIRST NODE OF LEVEL. UPDATE POINTERS. 828 : C 829 : LSTART = BACK + 1 830 : WIDTH = MAX0 (WIDTH, LSTART - LIST(NLEVEL)) 831 : IF ( WIDTH .GE. MAXWID ) GO TO 2000 832 : NLEVEL = NLEVEL - 1 833 : DEPTH = DEPTH + 1 834 : IF ( DEPTH .GT. MXDPTH ) GO TO 5000 835 : LIST (NLEVEL) = LSTART 836 : C 837 : C ... FIND ALL NEIGHBORS OF CURRENT NODE, ADD THEM TO QUEUE. 838 : C 839 : 1100 LFRONT = LIST (FRONT) 840 : SPTR = RSTART (LFRONT) 841 : FPTR = SPTR - DEGREE (LFRONT) - 1 842 : DO 1200 PTR = SPTR, FPTR 843 : NEWNOD = CONNEC (PTR) 844 : C 845 : C ... ADD TO QUEUE ONLY NODES NOT ALREADY IN QUEUE 846 : C 847 : IF ( DEGREE(NEWNOD) .LE. 0 ) GO TO 1200 848 : DEGREE (NEWNOD) = -DEGREE (NEWNOD) 849 : BACK = BACK + 1 850 : LIST (BACK) = NEWNOD 851 : 1200 CONTINUE 852 : FRONT = FRONT + 1 853 : C 854 : C ... IS QUEUE EMPTY? 855 : C ------------------- 856 : C 857 : IF ( FRONT .LE. BACK ) GO TO 1000 858 : C 859 : C ... YES, THE TREE IS BUILT. UNDO OUR MARKINGS. 860 : C 861 : IF (BACK .NE. ACTIVE) GO TO 6000 862 : C 863 : 1300 DO 1400 J = 1, BACK 864 : LISTJ = LIST(J) 865 : DEGREE (LISTJ) = -DEGREE (LISTJ) 866 : 1400 CONTINUE 867 : C 868 : RETURN 869 : C 870 : C ... ABORT GENERATION OF TREE BECAUSE IT IS ALREADY TOO WIDE 871 : C 872 : 2000 WIDTH = N + 1 873 : DEPTH = 0 874 : GO TO 1300 875 : C 876 : C ... INSUFFICIENT STORAGE ... 877 : C 878 : 5000 SPACE = 3 * ( (ACTIVE+1-BACK)*DEPTH / ACTIVE + (ACTIVE+1-BACK) ) 879 : ERROR = 111 880 : RETURN 881 : C 882 : 6000 ERROR = 13 883 : SPACE = -1 884 : RETURN 885 : C 886 : END 887 : SUBROUTINE GPSKCE (N, AVAIL, ACTIVE, DEPTH, WRKLEN, 888 : 1 LVLLST, LVLPTR, WORK, NXTNUM, TREE1, TREE2, 889 : 2 WIDTH1, WIDTH2, ONEIS1, ERROR, SPACE) 890 : C 891 : C ================================================================== 892 : C 893 : C TRANSITION BETWEEN ALGORITHM I AND ALGORITHM II OF 894 : C THE GIBBS-POOLE-STOCKMEYER PAPER. 895 : C 896 : C IN THIS IMPLEMENTATION ALGORITHM I REPRESENTS LEVEL TREES AS 897 : C LISTS OF NODES ORDERED BY LEVEL. ALGORITHM II APPEARS TO REQUIRE 898 : C LEVEL NUMBERS INDEXED BY NODE -- VECTORS FOR EFFICIENCY. 899 : C THIS SUBROUTINE CHANGES THE LEVEL TREE REPRESENTATION TO THAT 900 : C REQUIRED BY ALGORITHM II. NOTE THAT THE FIRST ALGORITHM CAN BE 901 : C CARRIED OUT WITH THE LEVEL NUMBER VECTOR FORMAT, PROBABLY REQURING 902 : C MORE COMPUTATION TIME, BUT PERHAPS LESS STORAGE. 903 : C 904 : C INPUT: TWO LEVEL TREES, AS LEVEL LISTS AND LEVEL POINTERS, 905 : C FOUND IN TWO OF THE THREE COLUMNS OF THE ARRAYS 'LVLLST' 906 : C AND 'LVLPTR' 907 : C 908 : C OUTPUT: TWO LEVEL TREES, AS VECTORS OF LEVEL NUMBERS, 909 : C ONE PACKED TO THE FRONT, ONE TO THE REAR OF THE WORKING 910 : C AREA 'WORK'. NOTE THAT 'WORK', 'LVLLST' AND 'LVLPTR' 911 : C SHARE COMMON LOCATIONS. 912 : C 913 : C ================================================================ 914 : C 915 : C ... STRUCTURE OF WORKSPACE 916 : C 917 : C INPUT .. (OUTPUT FROM GPSKCB) 918 : C 919 : C -------------------------------------------------------------- 920 : C : NUMBERED : TLIST1 PTR1 : TLIST2 PTR2 : TLIST3 PTR3 : 921 : C -------------------------------------------------------------- 922 : C 923 : C OUTPUT .. (GOES TO COMBIN) 924 : C 925 : C -------------------------------------------------------------- 926 : C : NUMBERED : TREE2 : ... : TREE1 : 927 : C -------------------------------------------------------------- 928 : C 929 : C ================================================================== 930 : C 931 : IMPLICIT INTEGER(I-N) 932 : INTEGER N, AVAIL, ACTIVE, DEPTH, WRKLEN, NXTNUM, 933 : 1 WIDTH1, WIDTH2, TREE1, TREE2, ERROR, SPACE 934 : C 935 : CIBM INTEGER *2 LVLLST(AVAIL,3), LVLPTR(AVAIL,3), WORK(WRKLEN) 936 : INTEGER LVLLST(AVAIL,3), LVLPTR(AVAIL,3), WORK(WRKLEN) 937 : C 938 : LOGICAL ONEIS1 939 : C 940 : C ------------------------------------------------------------------ 941 : C 942 : INTEGER I, BTREE, FTREE, FWIDTH, BWIDTH 943 : C 944 : C 945 : C ... CHECK THAT WE HAVE ENOUGH ROOM TO DO THE NECESSARY UNPACKING 946 : C 947 : IF (3*AVAIL .GT. WRKLEN) GO TO 6000 948 : IF (AVAIL .LT. N) GO TO 5100 949 : C 950 : C ... INPUT HAS THREE POSSIBLE CASES: 951 : C LVLLST(*,1) IS EMPTY 952 : C LVLLST(*,2) IS EMPTY 953 : C LVLLST(*,3) IS EMPTY 954 : C 955 : FTREE = TREE1 956 : BTREE = TREE2 957 : FWIDTH = WIDTH1 958 : BWIDTH = WIDTH2 959 : C 960 : TREE1 = WRKLEN - N + 1 961 : TREE2 = NXTNUM 962 : C 963 : IF ( (FTREE .EQ. 1) .OR. (BTREE .EQ. 1) ) GO TO 300 964 : C 965 : C ... CASE 1: 1ST SLOT IS EMPTY. UNPACK 3 INTO 1, 2 INTO 3 966 : C 967 : IF (FTREE .NE. 2) GO TO 100 968 : ONEIS1 = .TRUE. 969 : WIDTH2 = BWIDTH 970 : WIDTH1 = FWIDTH 971 : GO TO 200 972 : C 973 : 100 ONEIS1 = .FALSE. 974 : WIDTH1 = BWIDTH 975 : WIDTH2 = FWIDTH 976 : C 976 : 200 CALL GPSKCF (N, ACTIVE, DEPTH, LVLLST(1,3), LVLPTR(1,3), 977 : 1 WORK(TREE2), ONEIS1) 978 : C 978 : CALL GPSKCF (N, ACTIVE, DEPTH, LVLLST(1,2), LVLPTR(1,2), 979 : 1 WORK(TREE1), .NOT. ONEIS1) 980 : C 981 : GO TO 1000 982 : C 983 : C 984 : 300 IF ( (FTREE .EQ. 2) .OR. (BTREE .EQ. 2) ) GO TO 600 985 : C 986 : C ... CASE 2: 2ND SLOT IS EMPTY. TO ENABLE COMPLETE 987 : C REPACKING, MOVE 3 INTO 2, THEN FALL INTO NEXT CASE 988 : C 989 : DO 400 I = 1, ACTIVE 990 : LVLLST(I,2) = LVLLST(I,3) 991 : 400 CONTINUE 992 : C 993 : DO 500 I = 1, DEPTH 994 : LVLPTR(I,2) = LVLPTR(I,3) 995 : 500 CONTINUE 996 : C 997 : C ... CASE 3: SLOT 3 IS EMPTY. MOVE 1 INTO 3, THEN 2 INTO 1. 998 : C 999 : 600 IF (FTREE .EQ. 1) GO TO 700 1000 : ONEIS1 = .FALSE. 1001 : WIDTH1 = BWIDTH 1002 : WIDTH2 = FWIDTH 1003 : GO TO 800 1004 : C 1005 : 700 ONEIS1 = .TRUE. 1006 : WIDTH1 = FWIDTH 1007 : WIDTH2 = BWIDTH 1008 : C 1008 : 800 CALL GPSKCF (N, ACTIVE, DEPTH, LVLLST(1,1), LVLPTR(1,1), 1009 : 1 WORK(TREE1), .NOT. ONEIS1) 1010 : C 1010 : CALL GPSKCF (N, ACTIVE, DEPTH, LVLLST(1,2), LVLPTR(1,2), 1011 : 1 WORK(TREE2), ONEIS1) 1012 : 1000 RETURN 1013 : C 1014 : C ------------------------------------------------------------------ 1015 : C 1016 : 5100 SPACE = 3 * (N - AVAIL) 1017 : ERROR = 120 1018 : RETURN 1019 : C 1020 : 6000 ERROR = 20 1021 : SPACE = -1 1022 : RETURN 1023 : C 1024 : END 1025 : SUBROUTINE GPSKCF (N, ACTIVE, DEPTH, LVLLST, LVLPTR, LVLNUM, 1026 : 1 REVERS) 1027 : C 1028 : C ================================================================== 1029 : C 1030 : C CONVERT LEVEL STRUCTURE REPRESENTATION FROM A LIST OF NODES 1031 : C GROUPED BY LEVEL TO A VECTOR GIVING LEVEL NUMBER FOR EACH NODE. 1032 : C 1033 : C LVLLST, LVLPTR -- LIST OF LISTS 1034 : C 1035 : C LVLNUM -- OUTPUT VECTOR OF LEVEL NUMBERS 1036 : C 1037 : C REVERS -- IF .TRUE., NUMBER LEVEL STRUCTURE FROM BACK END 1038 : C INSTEAD OF FROM FRONT 1039 : C 1040 : C ================================================================== 1041 : C 1042 : IMPLICIT INTEGER(I-N) 1043 : INTEGER N, ACTIVE, DEPTH 1044 : C 1045 : CIBM INTEGER *2 LVLLST(ACTIVE), LVLPTR(DEPTH), LVLNUM(N) 1046 : INTEGER LVLLST(ACTIVE), LVLPTR(DEPTH), LVLNUM(N) 1047 : LOGICAL REVERS 1048 : C 1049 : C ------------------------------------------------------------------ 1050 : C 1051 : INTEGER I, LEVEL, LSTART, LEND, XLEVEL, PLSTRT, LVLLSI 1052 : C 1053 : IF (ACTIVE .EQ. N) GO TO 200 1054 : C 1055 : C ... IF NOT ALL NODES OF GRAPH ARE ACTIVE, MASK OUT THE 1056 : C NODES WHICH ARE NOT ACTIVE 1057 : C 1058 : DO 100 I = 1, N 1059 : LVLNUM(I) = 0 1060 : 100 CONTINUE 1061 : C 1062 : 200 DO 400 LEVEL = 1, DEPTH 1063 : XLEVEL = LEVEL 1064 : PLSTRT = DEPTH - LEVEL + 1 1065 : IF (REVERS) XLEVEL = PLSTRT 1066 : LSTART = LVLPTR (PLSTRT) 1067 : LEND = LVLPTR (PLSTRT - 1) - 1 1068 : C 1069 : DO 300 I = LSTART, LEND 1070 : LVLLSI = LVLLST(I) 1071 : LVLNUM (LVLLSI) = XLEVEL 1072 : 300 CONTINUE 1073 : 400 CONTINUE 1074 : C 1075 : RETURN 1076 : END 1077 : SUBROUTINE GPSKCG (N, DEGREE, RSTART, CONNEC, ACTIVE, WIDTH1, 1078 : 1 WIDTH2, TREE1, TREE2, WORK, WRKLEN, DEPTH, 1079 : 2 INC1, INC2, TOTAL, ONEIS1, REVRS1, ERROR, 1080 : 3 SPACE) 1081 : C 1082 : C ================================================================== 1083 : C 1084 : C COMBINE THE TWO ROOTED LEVEL TREES INTO A SINGLE LEVEL STRUCTURE 1085 : C WHICH MAY HAVE SMALLER WIDTH THAN EITHER OF THE TREES. THE NEW 1086 : C STRUCTURE IS NOT NECESSARILY A ROOTED STRUCTURE. 1087 : C 1088 : C PARAMETERS: 1089 : C 1090 : C N, DEGREE, RSTART, CONNEC -- GIVE THE DIMENSION AND STRUCTURE 1091 : C OF THE SPARSE SYMMETRIC MATRIX 1092 : C 1093 : C ACTIVE -- THE NUMBER OF NODES IN THIS CONNECTED COMPONENT OF 1094 : C THE MATRIX GRAPH 1095 : C 1096 : C TREE1 -- ON INPUT, ONE OF THE INPUT LEVEL TREES. ON 1097 : C OUTPUT, THE COMBINED LEVEL STRUCTURE 1098 : C 1099 : C TREE2 -- THE SECOND INPUT LEVEL TREE 1100 : C 1101 : C WIDTH1 -- THE MAXIMUM WIDTH OF A LEVEL IN TREE1 1102 : C 1103 : C WIDTH2 -- THE MAXIMUM WIDTH OF A LEVEL IN TREE2 1104 : C 1105 : C WORK -- A WORKING AREA OF LENGTH 'WRKLEN' 1106 : C 1107 : C INC1, -- VECTORS OF LENGTH 'DEPTH' 1108 : C INC2, 1109 : C TOTAL 1110 : C 1111 : C ONEIS1 -- INDICATES WHETHER TREE1 OR TREE2 REPRESENTS THE 1112 : C FORWARD TREE OR THE BACKWARDS TREE OF PHASE 1. 1113 : C USED TO MIMIC ARBITRARY TIE-BREAKING PROCEDURE OF 1114 : C ORIGINAL GIBBS-POOLE-STOCKMEYER CODE. 1115 : C 1116 : C REVRS1 -- OUTPUT PARAMETER INDICATING WHETHER A BACKWARDS 1117 : C ORDERING WAS USED FOR THE LARGEST COMPONENT OF 1118 : C THE REDUCED GRAPH 1119 : C 1120 : C ERROR -- NON-ZERO ONLY IF FAILURE OF SPACE ALLOCATION OR 1121 : C DATA STRUCTURE ERROR FOUND 1122 : C 1123 : C SPACE -- MINIMUM SPACE REQUIRED TO RERUN OR COMPLETE PHASE. 1124 : C 1125 : C ------------------------------------------------------------------ 1126 : C 1127 : IMPLICIT INTEGER(I-N) 1128 : INTEGER N, RSTART(N), ACTIVE, WIDTH1, WIDTH2, WRKLEN, DEPTH, 1129 : 2 ERROR, SPACE 1130 : C 1131 : CIBM INTEGER *2 DEGREE(N), CONNEC(1), TREE1(N), TREE2(N), 1132 : INTEGER DEGREE(N), CONNEC(1), TREE1(N), TREE2(N), 1133 : 1 WORK(WRKLEN), INC1(DEPTH), INC2(DEPTH), TOTAL(DEPTH) 1134 : C 1135 : LOGICAL ONEIS1, REVRS1 1136 : C 1137 : C ================================================================== 1138 : C 1139 : C << REMOVE ALL NODES OF PSEUDO-DIAMETERS >> 1140 : C << FIND CONNECTED COMPONENTS OF REDUCED GRAPH >> 1141 : C << COMBINE LEVEL TREES, COMPONENT BY COMPONENT >> 1142 : C 1143 : C ================================================================== 1144 : C 1145 : C STRUCTURE OF WORKSPACE ... 1146 : C 1147 : C ------------------------------------------------------------------ 1148 : C : NUMBERED : TREE2 : TOTAL : NODES : START : SIZE : INC1 : INC2 : 1149 : C ------------------------------------------------------------------ 1150 : C 1151 : C -------- 1152 : C TREE1 : 1153 : C -------- 1154 : C 1155 : C NUMBERED IS THE SET OF NUMBERED NODES (PROBABLY EMPTY) 1156 : C 1157 : C TREE1 AND TREE1 ARE LEVEL TREES (LENGTH N) 1158 : C TOTAL, INC1 AND INC2 ARE VECTORS OF NODE COUNTS PER LEVEL 1159 : C (LENGTH 'DEPTH') 1160 : C NODES IS THE SET OF NODES IN THE REDUCED GRAPH (THE NODES 1161 : C NOT ON ANY SHORTEST PATH FROM ONE END OF THE 1162 : C PSEUDODIAMETER TO THE OTHER) 1163 : C START, SIZE ARE POINTERS INTO 'NODES', ONE OF EACH FOR 1164 : C EACH CONNECTED COMPONENT OF THE REDUCED GRAPH. 1165 : C THE SIZES OF NODES, START AND SIZE ARE NOT KNOWN APRIORI. 1166 : C 1167 : C ================================================================== 1168 : INTEGER I, SIZE, AVAIL, CSTOP, START, COMPON, TREE1I, PCSTRT, 1169 : 1 CSTART, MXINC1, MXINC2, COMPNS, MXCOMP, OFFDIA, 1170 : 2 CSIZE, PCSIZE, WORKI, TWORKI 1171 : C 1172 : C ------------------------------------------------------------------ 1173 : C 1174 : C ... FIND ALL SHORTEST PATHS FROM START TO FINISH. REMOVE NODES ON 1175 : C THESE PATHS AND IN OTHER CONNECTED COMPONENTS OF FULL GRAPH 1176 : C FROM FURTHER CONSIDERATION. SIGN OF ENTRIES IN TREE1 IS USED 1177 : C AS A MASK. 1178 : C 1179 : OFFDIA = ACTIVE 1180 : C 1181 : DO 100 I = 1, DEPTH 1182 : TOTAL(I) = 0 1183 : 100 CONTINUE 1184 : C 1185 : DO 200 I = 1, N 1186 : TREE1I = TREE1 (I) 1187 : IF ((TREE1(I) .NE. TREE2(I)) .OR. (TREE1(I) .EQ. 0)) GO TO 200 1188 : TOTAL (TREE1I) = TOTAL (TREE1I) + 1 1189 : TREE1(I) = - TREE1(I) 1190 : OFFDIA = OFFDIA - 1 1191 : 200 CONTINUE 1192 : C 1193 : IF ( OFFDIA .EQ. 0 ) GO TO 1100 1194 : IF ( OFFDIA .LT. 0 ) GO TO 6000 1195 : C 1196 : C ... FIND CONNECTED COMPONENTS OF GRAPH INDUCED BY THE NODES NOT 1197 : C REMOVED. 'MXCOMP' IS THE LARGEST NUMBER OF COMPONENTS 1198 : C REPRESENTABLE IN THE WORKING SPACE AVAILABLE. 1199 : C 1200 : AVAIL = WRKLEN - OFFDIA 1201 : MXCOMP = AVAIL/2 1202 : START = OFFDIA + 1 1203 : SIZE = START + MXCOMP 1204 : C 1205 : IF (MXCOMP .LE. 0) GO TO 5100 1206 : C 1206 : CALL GPSKCH (N, DEGREE, RSTART, CONNEC, TREE1, OFFDIA, WORK, 1207 : 1 MXCOMP, WORK(START), WORK(SIZE), COMPNS, ERROR, 1208 : 2 SPACE) 1209 : IF ( ERROR .NE. 0 ) GO TO 5000 1210 : C 1211 : C ... RECORD SPACE ACTUALLY USED (NOT INCLUDING NUMBERED ) 1212 : C 1213 : SPACE = 2*N + 3*(DEPTH) + 2*COMPNS + OFFDIA 1214 : C 1215 : C ... SORT THE COMPONENT START POINTERS INTO INCREASING ORDER 1216 : C OF SIZE OF COMPONENT 1217 : C 1218 : IF (COMPNS .GT. 1) 1218 : 1 CALL GPSKCN (COMPNS, WORK(SIZE), WORK(START), ERROR) 1219 : IF (ERROR .NE. 0) GO TO 6200 1220 : C 1221 : C ... FOR EACH COMPONENT IN TURN, CHOOSE TO USE THE ORDERING OF THE 1222 : C 'FORWARD' TREE1 OR OF THE 'BACKWARD' TREE2 TO NUMBER THE NODES 1223 : C IN THIS COMPONENT. THE NUMBERING IS CHOSEN TO MINIMIZE THE 1224 : C MAXIMUM INCREMENT TO ANY LEVEL. 1225 : C 1226 : DO 1000 COMPON = 1, COMPNS 1227 : PCSTRT = START + COMPON - 1 1228 : CSTART = WORK (PCSTRT) 1229 : PCSIZE = SIZE + COMPON - 1 1230 : CSIZE = WORK (PCSIZE) 1231 : CSTOP = CSTART + CSIZE - 1 1232 : IF ( ( CSIZE .LT. 0 ) .OR. ( CSIZE .GT. OFFDIA ) ) GO TO 6100 1233 : C 1234 : DO 300 I = 1, DEPTH 1235 : INC1(I) = 0 1236 : INC2(I) = 0 1237 : 300 CONTINUE 1238 : C 1239 : MXINC1 = 0 1240 : MXINC2 = 0 1241 : C 1242 : DO 400 I = CSTART, CSTOP 1243 : WORKI = WORK(I) 1244 : TWORKI = -TREE1 (WORKI) 1245 : INC1 (TWORKI) = INC1 (TWORKI) + 1 1246 : TWORKI = TREE2 (WORKI) 1247 : INC2 (TWORKI) = INC2 (TWORKI) + 1 1248 : 400 CONTINUE 1249 : C 1250 : C ... BAROQUE TESTS BELOW DUPLICATE THE GIBBS-POOLE-STOCKMEYER- 1251 : C CRANE PROGRAM, *** NOT *** THE PUBLISHED ALGORITHM. 1252 : C 1253 : DO 500 I = 1, DEPTH 1254 : IF ((INC1(I) .EQ. 0) .AND. (INC2(I) .EQ. 0)) GO TO 500 1255 : IF (MXINC1 .LT. TOTAL(I) + INC1(I)) 1256 : 1 MXINC1 = TOTAL(I) + INC1(I) 1257 : IF (MXINC2 .LT. TOTAL(I) + INC2(I)) 1258 : 1 MXINC2 = TOTAL(I) + INC2(I) 1259 : 500 CONTINUE 1260 : C 1261 : C ... USE ORDERING OF NARROWER TREE UNLESS IT INCREASES 1262 : C WIDTH MORE THAN WIDER TREE. IN CASE OF TIE, USE TREE 2! 1263 : C 1264 : IF ( (MXINC1 .GT. MXINC2) .OR. 1265 : 1 ( (MXINC1 .EQ. MXINC2) .AND. ( (WIDTH1 .GT. WIDTH2) .OR. 1266 : 2 ( (WIDTH1 .EQ. WIDTH2) 1267 : 3 .AND. ONEIS1) ) ) ) 1268 : 4 GO TO 700 1269 : C 1270 : IF ( COMPON .EQ. 1 ) REVRS1 = .NOT. ONEIS1 1271 : C 1272 : DO 600 I = 1, DEPTH 1273 : TOTAL(I) = TOTAL(I) + INC1(I) 1274 : 600 CONTINUE 1275 : GO TO 1000 1276 : C 1277 : 700 IF ( COMPON .EQ. 1 ) REVRS1 = ONEIS1 1278 : DO 800 I = CSTART, CSTOP 1279 : WORKI = WORK(I) 1280 : TREE1 (WORKI) = - TREE2 (WORKI) 1281 : 800 CONTINUE 1282 : C 1283 : DO 900 I = 1, DEPTH 1284 : TOTAL(I) = TOTAL(I) + INC2(I) 1285 : 900 CONTINUE 1286 : C 1287 : 1000 CONTINUE 1288 : GO TO 2000 1289 : C 1290 : C ... DEFAULT WHEN THE REDUCED GRAPH IS EMPTY 1291 : C 1292 : 1100 REVRS1 = .TRUE. 1293 : SPACE = 2*N 1294 : C 1295 : 2000 RETURN 1296 : C 1297 : C ------------------------------------------------------------------ 1298 : C 1299 : C ERROR FOUND ... 1300 : C 1301 : 5000 SPACE = -1 1302 : GO TO 2000 1303 : C 1304 : 5100 SPACE = 2 - AVAIL 1305 : ERROR = 131 1306 : GO TO 2000 1307 : C 1308 : 6000 ERROR = 30 1309 : GO TO 5000 1310 : C 1311 : 6100 ERROR = 31 1312 : GO TO 5000 1313 : C 1314 : 6200 ERROR = 32 1315 : GO TO 5000 1316 : C 1317 : END 1318 : SUBROUTINE GPSKCH (N, DEGREE, RSTART, CONNEC, STATUS, NREDUC, 1319 : 1 WORK, MXCOMP, START, SIZE, COMPNS, ERROR, 1320 : 2 SPACE) 1321 : C 1322 : C ================================================================== 1323 : C 1324 : C FIND THE CONNECTED COMPONENTS OF THE GRAPH INDUCED BY THE SET 1325 : C OF NODES WITH POSITIVE 'STATUS'. WE SHALL BUILD THE LIST OF 1326 : C CONNECTED COMPONENTS IN 'WORK', WITH A LIST OF POINTERS 1327 : C TO THE BEGINNING NODES OF COMPONENTS LOCATED IN 'START' 1328 : C 1329 : C 1330 : IMPLICIT INTEGER(I-N) 1331 : INTEGER N, RSTART(N), NREDUC, MXCOMP, COMPNS, ERROR, SPACE 1332 : C 1333 : CIBM INTEGER *2 DEGREE(N), CONNEC(1), STATUS(N), WORK(NREDUC), 1334 : INTEGER DEGREE(N), CONNEC(1), STATUS(N), WORK(NREDUC), 1335 : 1 START(MXCOMP), SIZE(MXCOMP) 1336 : C 1337 : C 1338 : C PARAMETERS ... 1339 : C 1340 : C N -- DIMENSION OF THE ORIGINAL MATRIX 1341 : C DEGREE, RSTART, CONNEC -- THE STRUCTURE OF THE ORIGINAL MATRIX 1342 : C 1343 : C STATUS -- DERIVED FROM A LEVEL TREE. POSITIVE ENTRIES INDICATE 1344 : C ACTIVE NODES. NODES WITH STATUS <= 0 ARE IGNORED. 1345 : C 1346 : C NREDUC -- THE NUMBER OF ACTIVE NODES 1347 : C 1348 : C WORK -- WORK SPACE, USED AS A QUEUE TO BUILD CONNECTED 1349 : C COMPONENTS IN PLACE. 1350 : C 1351 : C MXCOMP -- MAXIMUM NUMBER OF COMPONENTS ALLOWED BY CURRENT 1352 : C SPACE ALLOCATION. MUST NOT BE VIOLATED. 1353 : C 1354 : C START -- POINTER TO BEGINNING OF I-TH CONNECTED COMPONENT 1355 : C 1356 : C SIZE -- SIZE OF EACH COMPONENT 1357 : C 1358 : C COMPNS -- NUMBER OF COMPONENTS ACTUALLY FOUND 1359 : C 1360 : C ERROR -- SHOULD BE ZERO ON RETURN UNLESS WE HAVE TOO LITTLE 1361 : C SPACE OR WE ENCOUNTER AN ERROR IN THE DATA STRUCTURE 1362 : C 1363 : C SPACE -- MAXIMUM AMOUNT OF WORKSPACE USED / NEEDED 1364 : C 1365 : C ================================================================== 1366 : C 1367 : INTEGER I, J, FREE, JPTR, NODE, JNODE, FRONT, CDGREE, ROOT 1368 : C 1369 : C ------------------------------------------------------------------ 1370 : C 1371 : C 1372 : C REPEAT 1373 : C << FIND AN UNASSIGNED NODE AND START A NEW COMPONENT >> 1374 : C REPEAT 1375 : C << ADD ALL NEW NEIGHBORS OF FRONT NODE TO QUEUE, >> 1376 : C << REMOVE FRONT NODE. >> 1377 : C UNTIL <> 1378 : C UNTIL << ALL NODES ASSIGNED >> 1379 : C 1380 : FREE = 1 1381 : COMPNS = 0 1382 : ROOT = 1 1383 : C 1384 : C ... START OF OUTER REPEAT LOOP 1385 : C 1386 : C ... FIND AN UNASSIGNED NODE 1387 : C 1388 : 100 DO 200 I = ROOT, N 1389 : IF (STATUS(I) .LE. 0) GO TO 200 1390 : NODE = I 1391 : GO TO 300 1392 : 200 CONTINUE 1393 : GO TO 6100 1394 : C 1395 : C ... START NEW COMPONENT 1396 : C 1397 : 300 COMPNS = COMPNS + 1 1398 : ROOT = NODE + 1 1399 : IF (COMPNS .GT. MXCOMP) GO TO 5000 1400 : START (COMPNS) = FREE 1401 : WORK (FREE) = NODE 1402 : STATUS (NODE) = -STATUS (NODE) 1403 : FRONT = FREE 1404 : FREE = FREE + 1 1405 : C 1406 : C ... INNER REPEAT UNTIL QUEUE BECOMES EMPTY 1407 : C 1408 : 400 NODE = WORK (FRONT) 1409 : FRONT = FRONT + 1 1410 : C 1411 : JPTR = RSTART (NODE) 1412 : CDGREE = DEGREE (NODE) 1413 : DO 500 J = 1, CDGREE 1414 : JNODE = CONNEC (JPTR) 1415 : JPTR = JPTR + 1 1416 : IF (STATUS(JNODE) .LT. 0) GO TO 500 1417 : IF (STATUS(JNODE) .EQ. 0) GO TO 6000 1418 : STATUS (JNODE) = -STATUS (JNODE) 1419 : WORK (FREE) = JNODE 1420 : FREE = FREE + 1 1421 : 500 CONTINUE 1422 : C 1423 : IF (FRONT .LT. FREE) GO TO 400 1424 : C 1425 : C ... END OF INNER REPEAT. COMPUTE SIZE OF COMPONENT AND 1426 : C SEE IF THERE ARE MORE NODES TO BE ASSIGNED 1427 : C 1428 : SIZE (COMPNS) = FREE - START (COMPNS) 1429 : IF (FREE .LE. NREDUC) GO TO 100 1430 : C 1431 : IF (FREE .NE. NREDUC+1) GO TO 6200 1432 : RETURN 1433 : C 1434 : C ------------------------------------------------------------------ 1435 : C 1436 : 5000 SPACE = NREDUC - FREE + 1 1437 : ERROR = 130 1438 : RETURN 1439 : C 1440 : 6000 ERROR = 33 1441 : SPACE = -1 1442 : RETURN 1443 : C 1444 : 6100 ERROR = 34 1445 : SPACE = -1 1446 : RETURN 1447 : C 1448 : 6200 ERROR = 35 1449 : SPACE = -1 1450 : RETURN 1451 : END 1452 : SUBROUTINE GPSKCI (N, ACTIVE, DEPTH, LSTRUC, LVLLST, LVLPTR, 1453 : 1 LTOTAL, ERROR, SPACE) 1454 : C 1455 : C ================================================================== 1456 : C 1457 : C TRANSITIONAL SUBROUTINE, ALGORITHM II TO IIIA OR IIIB. 1458 : C 1459 : C CONVERT LEVEL STRUCTURE GIVEN AS VECTOR OF LEVEL NUMBERS FOR NODES 1460 : C TO STRUCTURE AS LIST OF NODES BY LEVEL 1461 : C 1462 : C N, ACTIVE, DEPTH -- PROBLEM SIZES 1463 : C LSTRUC -- INPUT LEVEL STRUCTURE 1464 : C LVLLST, LVLPTR -- OUTPUT LEVEL STRUCTURE 1465 : C LTOTAL -- NUMBER OF NODES AT EACH LEVEL (PRECOMPUTED) 1466 : C 1467 : IMPLICIT INTEGER(I-N) 1468 : INTEGER N, ACTIVE, DEPTH, ERROR, SPACE 1469 : C 1470 : CIBM INTEGER *2 LSTRUC(N), LVLLST(ACTIVE), LVLPTR(1), LTOTAL(DEPTH) 1471 : INTEGER LSTRUC(N), LVLLST(ACTIVE), LVLPTR(1), LTOTAL(DEPTH) 1472 : C 1473 : C =============================================================== 1474 : C 1475 : C STRUCTURE OF WORKSPACE .. 1476 : C 1477 : C INPUT (FROM COMBIN) .. 1478 : C 1479 : C ------------------------------------------------------------------ 1480 : C : NUMBERED : ..(N).. : TOTAL : ... : TREE : 1481 : C ------------------------------------------------------------------ 1482 : C 1483 : C OUTPUT (TO GPSKCJ OR GPSKCK) .. 1484 : C 1485 : C ------------------------------------------------------------------ 1486 : C : NUMBERED : ... : TLIST : TPTR : TREE : 1487 : C ------------------------------------------------------------------ 1488 : C 1489 : C HERE, NUMBERED IS THE SET OF NODES IN NUMBERED COMPONENTS 1490 : C TOTAL IS A VECTOR OF LENGTH 'DEPTH' GIVING THE NUMBER 1491 : C OF NODES IN EACH LEVEL OF THE 'TREE'. 1492 : C TLIST, TPTR ARE LISTS OF NODES OF THE TREE, ARRANGED 1493 : C BY LEVEL. TLIST IS OF LENGTH 'ACTIVE', TPTR 'DEPTH+1'. 1494 : C 1495 : C ================================================================= 1496 : C 1497 : INTEGER I, ACOUNT, START, LEVEL, PLEVEL 1498 : C 1499 : C ... ESTABLISH STARTING AND ENDING POINTERS FOR EACH LEVEL 1500 : C 1501 : START = 1 1502 : DO 100 I = 1, DEPTH 1503 : LVLPTR(I) = START 1504 : START = START + LTOTAL(I) 1505 : LTOTAL(I) = START 1506 : 100 CONTINUE 1507 : LVLPTR(DEPTH+1) = START 1508 : C 1509 : ACOUNT = 0 1510 : DO 300 I = 1, N 1511 : IF (LSTRUC(I)) 200, 300, 6000 1512 : 200 LEVEL = -LSTRUC(I) 1513 : LSTRUC(I) = LEVEL 1514 : PLEVEL = LVLPTR (LEVEL) 1515 : LVLLST (PLEVEL) = I 1516 : LVLPTR (LEVEL) = LVLPTR (LEVEL) + 1 1517 : ACOUNT = ACOUNT + 1 1518 : IF (LVLPTR (LEVEL) .GT. LTOTAL (LEVEL)) GO TO 6100 1519 : 300 CONTINUE 1520 : C 1521 : C ... RESET STARTING POINTERS 1522 : C 1523 : LVLPTR(1) = 1 1524 : DO 400 I = 1, DEPTH 1525 : LVLPTR(I+1) = LTOTAL(I) 1526 : 400 CONTINUE 1527 : C 1528 : RETURN 1529 : C 1530 : C ------------------------------------------------------------------ 1531 : C 1532 : 6000 ERROR = 40 1533 : GO TO 6200 1534 : C 1535 : 6100 ERROR = 41 1536 : C 1537 : 6200 SPACE = -1 1538 : RETURN 1539 : C 1540 : END 1541 : SUBROUTINE GPSKCJ (N, DEGREE, RSTART, CONNEC, 1542 : 1 NCOMPN, INVNUM, SNODE1, SNODE2, REVRS1, 1543 : 2 DEPTH, LVLLST, LVLPTR, LVLNUM, ERROR, 1544 : 3 SPACE) 1545 : C 1546 : C ================================================================== 1547 : C 1548 : C NUMBER THE NODES IN A GENERALIZED LEVEL STRUCTURE ACCORDING 1549 : C TO A GENERALIZATION OF THE CUTHILL MCKEE STRATEGY. 1550 : C 1551 : C N -- DIMENSION OF ORIGINAL PROBLEM 1552 : C DEGREE, RSTART, CONNEC -- GIVE STRUCTURE OF SPARSE AND 1553 : C SYMMETRIC MATRIX 1554 : C 1555 : C NCOMPN -- NUMBER OF NODES IN THIS COMPONENT OF MATRIX GRAPH 1556 : C 1557 : C INVNUM -- WILL BECOME A LIST OF THE ORIGINAL NODES IN THE ORDER 1558 : C WHICH REDUCES THE BANDWIDTH OF THE MATRIX. 1559 : C 1560 : C NXTNUM -- THE NEXT INDEX TO BE ASSIGNED (1 FOR FIRST COMPONENT) 1561 : C 1562 : C REVRS1 -- IF .TRUE., FIRST COMPONENT OF REDUCED GRAPH WAS NUMBERED 1563 : C BACKWARDS. 1564 : C 1565 : C LVLLST -- LIST OF NODES IN LEVEL TREE ORDERED BY LEVEL. 1566 : C 1567 : C LVLPTR -- POSITION OF INITIAL NODE IN EACH LEVEL OF LVLLST. 1568 : C 1569 : C LVLNUM -- LEVEL NUMBER OF EACH NODE IN COMPONENT 1570 : C 1571 : C 1572 : IMPLICIT INTEGER(I-N) 1573 : INTEGER N, RSTART(N), NCOMPN, SNODE1, SNODE2, DEPTH, 1574 : 1 ERROR, SPACE 1575 : C 1576 : CIBM INTEGER *2 DEGREE(N), CONNEC(1), INVNUM(NCOMPN), 1577 : INTEGER DEGREE(N), CONNEC(1), INVNUM(NCOMPN), 1578 : 1 LVLLST(NCOMPN), LVLPTR(DEPTH), LVLNUM(N) 1579 : C 1580 : LOGICAL REVRS1 1581 : C 1582 : C 1583 : C ================================================================== 1584 : C 1585 : C NUMBERING REQUIRES TWO QUEUES, WHICH CAN BE BUILD IN PLACE 1586 : C IN INVNUM. 1587 : C 1588 : C 1589 : C ================================================================== 1590 : C A L G O R I T H M S T R U C T U R E 1591 : C ================================================================== 1592 : C 1593 : C << SET QUEUE1 TO BE THE SET CONTAINING ONLY THE START NODE. >> 1594 : C 1595 : C FOR LEVEL = 1 TO DEPTH DO 1596 : C 1597 : C BEGIN 1598 : C LOOP 1599 : C 1600 : C REPEAT 1601 : C BEGIN 1602 : C << CNODE <- FRONT OF QUEUE1 >> 1603 : C << ADD UNNUMBERED NEIGHBORS OF CNODE TO THE BACK >> 1604 : C << OF QUEUE1 OR QUEUE2 (USE QUEUE1 IF NEIGHBOR >> 1605 : C << AT SAME LEVEL, QUEUE2 IF AT NEXT LEVEL). SORT >> 1606 : C << THE NEWLY QUEUED NODES INTO INCREASING ORDER OF >> 1607 : C << DEGREE. NUMBER CNODE, DELETE IT FROM QUEUE1. >> 1608 : C END 1609 : C UNTIL 1610 : C << QUEUE1 IS EMPTY >> 1611 : C 1612 : C EXIT IF << ALL NODES AT THIS LEVEL NUMBERED >> 1613 : C 1614 : C BEGIN 1615 : C << FIND THE UNNUMBERED NODE OF MINIMAL DEGREE AT THIS >> 1616 : C << LEVEL, RESTART QUEUE1 WITH THIS NODE. >> 1617 : C END 1618 : C 1619 : C END << LOOP LOOP >> 1620 : C 1621 : C << PROMOTE QUEUE2 TO BE INITIAL QUEUE1 FOR NEXT ITERATION >> 1622 : C << OF FOR LOOP. >> 1623 : C 1624 : C END < > 1625 : C 1626 : C ================================================================== 1627 : C 1628 : C STRUCTURE OF WORKSPACE .. 1629 : C 1630 : C -------------------------------------------------------------- 1631 : C : NUMBERED : QUEUE1 : QUEUE2 : ... : TLIST : TPTR : TREE : 1632 : C -------------------------------------------------------------- 1633 : C 1634 : C ON COMPLETION, WE HAVE ONLY A NEW, LONGER NUMBERED SET. 1635 : C 1636 : C ================================================================== 1637 : INTEGER I, BQ1, BQ2, FQ1, INC, CPTR, CNODE, 1638 : 1 INODE, LEVEL, NLEFT, LSTART, LWIDTH, QUEUE1, 1639 : 2 QUEUE2, CDGREE, XLEVEL, STNODE, ILEVEL, SQ1, SQ2, 1640 : 3 NSORT, LOWDG, BPTR, LVLLSC, LVLLSB, INVNMI 1641 : C 1642 : LOGICAL FORWRD, RLEVEL 1643 : C 1644 : C ------------------------------------------------------------------ 1645 : C 1646 : C ... GIBBS-POOLE-STOCKMEYER HEURISTIC CHOICE OF ORDER 1647 : C 1648 : IF (DEGREE(SNODE1) .GT. DEGREE(SNODE2)) GO TO 10 1649 : FORWRD = REVRS1 1650 : STNODE = SNODE1 1651 : GO TO 20 1652 : C 1653 : 10 FORWRD = .NOT. REVRS1 1654 : STNODE = SNODE2 1655 : C 1656 : C ... SET UP INITIAL QUEUES AT FRONT OF 'INVNUM' FOR FORWRD ORDER, 1657 : C AT BACK FOR REVERSED ORDER. 1658 : C 1659 : 20 IF (FORWRD) GO TO 100 1660 : INC = -1 1661 : QUEUE1 = NCOMPN 1662 : GO TO 200 1663 : C 1664 : 100 INC = +1 1665 : QUEUE1 = 1 1666 : C 1667 : 200 INVNUM (QUEUE1) = STNODE 1668 : RLEVEL = (LVLNUM(STNODE) .EQ. DEPTH) 1669 : LVLNUM (STNODE) = 0 1670 : FQ1 = QUEUE1 1671 : BQ1 = QUEUE1 + INC 1672 : C 1673 : C ------------------------------- 1674 : C NUMBER NODES LEVEL BY LEVEL ... 1675 : C ------------------------------- 1676 : C 1677 : DO 3000 XLEVEL = 1, DEPTH 1678 : LEVEL = XLEVEL 1679 : IF (RLEVEL) LEVEL = DEPTH - XLEVEL + 1 1680 : C 1681 : LSTART = LVLPTR (LEVEL) 1682 : LWIDTH = LVLPTR (LEVEL+1) - LSTART 1683 : NLEFT = LWIDTH 1684 : QUEUE2 = QUEUE1 + INC*LWIDTH 1685 : BQ2 = QUEUE2 1686 : C 1687 : C ============================================================== 1688 : C ... 'LOOP' CONSTRUCT BEGINS AT STATEMENT 1000 1689 : C THE INNER 'REPEAT' WILL BE DONE AS MANY TIMES AS 1690 : C IS NECESSARY TO NUMBER ALL THE NODES AT THIS LEVEL. 1691 : C ============================================================== 1692 : C 1693 : 1000 CONTINUE 1694 : C 1695 : C ========================================================== 1696 : C ... REPEAT ... UNTIL QUEUE1 BECOMES EMPTY 1697 : C TAKE NODE FROM FRONT OF QUEUE1, FIND EACH OF ITS 1698 : C NEIGHBORS WHICH HAVE NOT YET BEEN NUMBERED, AND 1699 : C ADD THE NEIGHBORS TO QUEUE1 OR QUEUE2 ACCORDING TO 1700 : C THEIR LEVELS. 1701 : C ========================================================== 1702 : C 1703 : 1100 CNODE = INVNUM (FQ1) 1704 : FQ1 = FQ1 + INC 1705 : SQ1 = BQ1 1706 : SQ2 = BQ2 1707 : NLEFT = NLEFT - 1 1708 : C 1709 : CPTR = RSTART (CNODE) 1710 : CDGREE = DEGREE (CNODE) 1711 : DO 1300 I = 1, CDGREE 1712 : INODE = CONNEC (CPTR) 1713 : CPTR = CPTR + 1 1714 : ILEVEL = LVLNUM (INODE) 1715 : IF (ILEVEL .EQ. 0) GO TO 1300 1716 : LVLNUM (INODE) = 0 1717 : IF ( ILEVEL .EQ. LEVEL ) GO TO 1200 1718 : C 1719 : IF (IABS(LEVEL-ILEVEL) .NE. 1) GO TO 6400 1720 : INVNUM (BQ2) = INODE 1721 : BQ2 = BQ2 + INC 1722 : GO TO 1300 1723 : C 1724 : 1200 INVNUM (BQ1) = INODE 1725 : BQ1 = BQ1 + INC 1726 : 1300 CONTINUE 1727 : C 1728 : C ================================================== 1729 : C ... SORT THE NODES JUST ADDED TO QUEUE1 AND QUEUE2 1730 : C SEPARATELY INTO INCREASING ORDER OF DEGREE. 1731 : C ================================================== 1732 : C 1733 : IF (IABS (BQ1 - SQ1) .LE. 1) GO TO 1500 1734 : NSORT = IABS (BQ1 - SQ1) 1735 : IF (FORWRD) GO TO 1400 1735 : CALL GPSKCP (NSORT, INVNUM(BQ1+1), N, DEGREE, 1736 : 1 ERROR) 1737 : IF (ERROR .NE. 0) GO TO 6600 1738 : GO TO 1500 1739 : C 1739 : 1400 CALL GPSKCQ (NSORT, INVNUM(SQ1), N, DEGREE, 1740 : 1 ERROR) 1741 : IF (ERROR .NE. 0) GO TO 6600 1742 : C 1743 : 1500 IF (IABS (BQ2 - SQ2) .LE. 1) GO TO 1700 1744 : NSORT = IABS (BQ2 - SQ2) 1745 : IF (FORWRD) GO TO 1600 1745 : CALL GPSKCP (NSORT, INVNUM(BQ2+1), N, DEGREE, 1746 : 1 ERROR) 1747 : IF (ERROR .NE. 0) GO TO 6600 1748 : GO TO 1700 1749 : C 1749 : 1600 CALL GPSKCQ (NSORT, INVNUM(SQ2), N, DEGREE, 1750 : 1 ERROR) 1751 : IF (ERROR .NE. 0) GO TO 6600 1752 : C 1753 : C ... END OF REPEAT LOOP 1754 : C 1755 : 1700 IF (FQ1 .NE. BQ1) GO TO 1100 1756 : C 1757 : C ============================================================== 1758 : C ... QUEUE1 IS NOW EMPTY ... 1759 : C IF THERE ARE ANY UNNUMBERED NODES LEFT AT THIS LEVEL, 1760 : C FIND THE ONE OF MINIMAL DEGREE AND RETURN TO THE 1761 : C REPEAT LOOP ABOVE. 1762 : C ============================================================== 1763 : C 1764 : 2000 IF ((BQ1 .EQ. QUEUE2) .AND. (NLEFT .EQ. 0)) GO TO 2900 1765 : C 1766 : IF ((NLEFT .LE. 0) .OR. (NLEFT .NE. INC * (QUEUE2 - BQ1))) 1767 : 1 GO TO 6200 1768 : C 1769 : LOWDG = N + 1 1770 : BPTR = N + 1 1771 : CPTR = LSTART - 1 1772 : DO 2800 I = 1, NLEFT 1773 : 2600 CPTR = CPTR + 1 1774 : LVLLSC = LVLLST (CPTR) 1775 : IF (LVLNUM (LVLLSC) .EQ. LEVEL) GO TO 2700 1776 : IF (LVLNUM (LVLLSC) .NE. 0) GO TO 6300 1777 : GO TO 2600 1778 : C 1779 : 2700 IF (DEGREE(LVLLSC) .GE. LOWDG) GO TO 2800 1780 : LOWDG = DEGREE (LVLLSC) 1781 : BPTR = CPTR 1782 : C 1783 : 2800 CONTINUE 1784 : C 1785 : C ... MINIMAL DEGREE UNNUMBERED NODE FOUND ... 1786 : C 1787 : IF (BPTR .GT. N) GO TO 6500 1788 : LVLLSB = LVLLST (BPTR) 1789 : INVNUM (BQ1) = LVLLSB 1790 : LVLNUM (LVLLSB) = 0 1791 : BQ1 = BQ1 + INC 1792 : GO TO 1000 1793 : C 1794 : C ============================================= 1795 : C ... ADVANCE QUEUE POINTERS TO MAKE QUEUE2 THE 1796 : C NEW QUEUE1 FOR THE NEXT ITERATION. 1797 : C ============================================= 1798 : C 1799 : 2900 QUEUE1 = QUEUE2 1800 : FQ1 = QUEUE1 1801 : BQ1 = BQ2 1802 : IF ((BQ1 .EQ. FQ1) .AND. (XLEVEL .LT. DEPTH)) GO TO 6100 1803 : C 1804 : 3000 CONTINUE 1805 : C 1806 : C ... CHANGE SIGN OF DEGREE TO MARK THESE NODES AS 'NUMBERED' 1807 : C 1808 : DO 3100 I = 1, NCOMPN 1809 : INVNMI = INVNUM(I) 1810 : DEGREE (INVNMI) = -DEGREE (INVNMI) 1811 : 3100 CONTINUE 1812 : C 1813 : RETURN 1814 : C 1815 : C ------------------------------------------------------------------ 1816 : C 1817 : 6000 SPACE = -1 1818 : RETURN 1819 : C 1820 : 6100 ERROR = 51 1821 : GO TO 6000 1822 : C 1823 : 6200 ERROR = 52 1824 : GO TO 6000 1825 : C 1826 : 6300 ERROR = 53 1827 : GO TO 6000 1828 : C 1829 : 6400 ERROR = 54 1830 : GO TO 6000 1831 : C 1832 : 6500 ERROR = 55 1833 : GO TO 6000 1834 : C 1835 : 6600 ERROR = 56 1836 : GO TO 6000 1837 : C 1838 : END 1839 : SUBROUTINE GPSKCK (N, DEGREE, RSTART, CONNEC, WRKLEN, NXTNUM, 1840 : 1 WORK, NCOMPN, DEPTH, LVLLST, LVLPTR, LVLNUM, 1841 : 2 ERROR, SPACE) 1842 : C 1843 : IMPLICIT INTEGER(I-N) 1844 : INTEGER N, RSTART(N), WRKLEN, NXTNUM, NCOMPN, DEPTH, ERROR, 1845 : 1 SPACE 1846 : C 1847 : CIBM INTEGER *2 DEGREE(N), CONNEC(1), WORK(WRKLEN), LVLLST(N), 1848 : INTEGER DEGREE(N), CONNEC(1), WORK(WRKLEN), LVLLST(N), 1849 : 1 LVLPTR(DEPTH), LVLNUM(N) 1850 : C 1851 : C ================================================================== 1852 : C 1853 : C NUMBER NODES IN A GENERALIZED LEVEL STRUCTURE ACCORDING TO 1854 : C A GENERALIZATION OF THE KING ALGORITHM, WHICH REDUCES 1855 : C THE PROFILE OF THE SPARSE SYMMETRIC MATRIX. 1856 : C 1857 : C --------------------- 1858 : C 1859 : C CODE USES A PRIORITY QUEUE TO CHOOSE THE NEXT NODE TO BE NUMBERED 1860 : C THE PRIORITY QUEUE IS REPRESENTED BY A SIMPLE LINEAR-LINKED LIST 1861 : C TO SAVE SPACE. THIS WILL REQUIRE MORE SEARCHING THAN A FULLY 1862 : C LINKED REPRESENTATION, BUT THE DATA MANIPULATION IS SIMPLER. 1863 : C 1864 : C ------------------- 1865 : C 1866 : C << ESTABLISH PRIORITY QUEUE 'ACTIVE' FOR LEVEL 1 NODES >> 1867 : C 1868 : C FOR I = 1 TO DEPTH DO 1869 : C << SET QUEUE 'QUEUED' TO BE EMPTY, LIST 'NEXT' TO BE >> 1870 : C << SET OF NODES AT NEXT LEVEL. >> 1871 : C 1872 : C FOR J = 1 TO 'NODES AT THIS LEVEL' DO 1873 : C << FIND FIRST NODE IN ACTIVE WITH MINIMAL CONNECTIONS >> 1874 : C << TO 'NEXT'. NUMBER THIS NODE AND REMOVE HIM FROM >> 1875 : C << 'ACTIVE'. FOR EACH NODE IN 'NEXT' WHICH CONNECTED >> 1876 : C << TO THIS NODE, MOVE IT TO 'QUEUED' AND REMOVE IT >> 1877 : C << FROM 'NEXT'. >> 1878 : C 1879 : C << SET NEW QUEUE 'ACTIVE' TO BE 'QUEUED' FOLLOWED BY ANY >> 1880 : C << NODES STILL IN 'NEXT'. >> 1881 : C 1882 : C ================================================================== 1883 : C 1884 : C DATA STRUCTURE ASSUMPTIONS: 1885 : C THE FIRST 'NXTNUM-1' ELEMENTS OF WORK ARE ALREADY IN USE. 1886 : C THE LEVEL STRUCTURE 'LVLLST' IS CONTIGUOUS WITH WORK, THAT IS, 1887 : C IT RESIDES IN ELEMENTS WRKLEN+1, ... OF WORK. 'LVLPTR' AND 1888 : C 'LVLNUM' ARE ALSO EMBEDDED IN WORK, BEHIND 'LVLLST'. THE 1889 : C THREE VECTORS ARE PASSED SEPARATELY TO CLARIFY THE INDEXING, 1890 : C BUT THE QUEUES DEVELOPED WILL BE ALLOWED TO OVERRUN 'LVLLST' 1891 : C AS NEEDED. 1892 : C 1893 : C ... BUILD THE FIRST 'ACTIVE' QUEUE STARTING W1 LOCATIONS FROM 1894 : C THE FRONT OF THE CURRENT WORKING AREA (W1 IS THE WIDTH OF THE 1895 : C FIRST LEVEL). BUILD THE FIRST 'QUEUED' QUEUE STARTING FROM 1896 : C THE BACK OF WORK SPACE. THE LIST 'NEXT' WILL BE REALIZED 1897 : C LVLNUM(I) > 0 <== LEVEL NUMBER OF NODE. 'NEXT' IS 1898 : C SET WITH LVLNUM(I) = LEVEL+1 1899 : C LVLNUM(I) = 0 <== I-TH NODE IS IN 'QUEUED' OR IS 1900 : C NOT IN THIS COMPONENT OF GRAPH, 1901 : C OR HAS JUST BEEN NUMBERED. 1902 : C LVLNUM(I) < 0 <== I-TH NODE IS IN 'ACTIVE' AND IS 1903 : C CONNECTED TO -LVLNUM(I) NODES IN 1904 : C 'NEXT'. 1905 : C 1906 : C ================================================================== 1907 : C 1908 : C STRUCTURE OF WORKSPACE .. 1909 : C 1910 : C -------------------------------------------------------------- 1911 : C : NUMBERED : DONE : ACTIVE : ALEVEL : ... : QUEUED : LVLLST : 1912 : C -------------------------------------------------------------- 1913 : C 1914 : C ------------------- 1915 : C LVLPTR : LVLNUM : 1916 : C ------------------- 1917 : C 1918 : C IN THE ABOVE, 1919 : C NUMBERED IS THE SET OF NODES ALREADY NUMBERED FROM 1920 : C PREVIOUS COMPONENTS AND EARLIER LEVELS OF THIS COMPONENT. 1921 : C DONE, ACTIVE, ALEVEL ARE VECTORS OF LENGTH THE WIDTH OF 1922 : C THE CURRENT LEVEL. ACTIVE IS A SET OF INDICES INTO 1923 : C ALEVEL. AS THE NODES IN ALEVEL ARE NUMBERED, THEY 1924 : C ARE PLACED INTO 'DONE'. 1925 : C QUEUED IS A QUEUE OF NODES IN THE 'NEXT' LEVEL, WHICH 1926 : C GROWS FROM THE START OF THE 'NEXT' LEVEL IN LVLLST 1927 : C FORWARDS TOWARD 'ALEVEL'. QUEUED IS OF LENGTH NO MORE 1928 : C THAN THE WIDTH OF THE NEXT LEVEL. 1929 : C LVLLST IS THE LIST OF UNNUMBERED NODES IN THE TREE, 1930 : C ARRANGED BY LEVEL. 1931 : C 1932 : C ================================================================== 1933 : INTEGER I, J, K, PTR, JPTR, KPTR, LPTR, MPTR, PPTR, RPTR, 1934 : 1 MPPTR, JNODE, KNODE, CNODE, LEVEL, LOWDG, UNUSED, 1935 : 2 MXQUE, NNEXT, ASTART, MINDG, LSTART, LWIDTH, ACTIVE, 1936 : 2 QUEUEB, QUEUED, QCOUNT, NCONNC, NACTIV, CDGREE, 1937 : 3 LDGREE, NFINAL, JDGREE, STRTIC, ADDED, TWRKLN, 1938 : 4 LVLLSL, CONNEJ, CONNER, ASTPTR, ACTPTR, ACTIVI, 1939 : 5 ASTRTI, QUEUEI, ACPPTR 1940 : C 1941 : C ------------------------------------------------------------------ 1942 : C 1943 : TWRKLN = WRKLEN + NCOMPN + N + DEPTH + 1 1944 : UNUSED = TWRKLN 1945 : C 1946 : ASTART = LVLPTR(1) 1947 : LWIDTH = LVLPTR(2) - ASTART 1948 : ASTART = WRKLEN + 1 1949 : ACTIVE = NXTNUM + LWIDTH + 1 1950 : NACTIV = LWIDTH 1951 : NFINAL = NXTNUM + NCOMPN 1952 : C 1953 : NNEXT = LVLPTR(3) - LVLPTR(2) 1954 : QUEUED = WRKLEN 1955 : QUEUEB = QUEUED 1956 : MXQUE = ACTIVE + LWIDTH 1957 : C 1958 : C ... BUILD FIRST PRIORITY QUEUE 'ACTIVE' 1959 : C 1960 : LOWDG = - (N + 1) 1961 : LPTR = LVLPTR(1) 1962 : DO 200 I = 1, LWIDTH 1963 : NCONNC = 0 1964 : LVLLSL= LVLLST (LPTR) 1965 : JPTR = RSTART (LVLLSL) 1966 : LDGREE = DEGREE(LVLLSL) 1967 : DO 100 J = 1, LDGREE 1968 : CONNEJ = CONNEC (JPTR) 1969 : IF ( LVLNUM (CONNEJ) .EQ. 2 ) NCONNC = NCONNC - 1 1970 : JPTR = JPTR + 1 1971 : 100 CONTINUE 1972 : C 1973 : ACTIVI = ACTIVE + I - 1 1974 : WORK (ACTIVI) = I 1975 : LVLNUM (LVLLSL) = NCONNC 1976 : LOWDG = MAX0 (LOWDG, NCONNC) 1977 : LPTR = LPTR + 1 1978 : 200 CONTINUE 1979 : WORK (ACTIVE-1) = 0 1980 : C 1981 : C ----------------------------------- 1982 : C NOW NUMBER NODES LEVEL BY LEVEL ... 1983 : C ----------------------------------- 1984 : C 1985 : DO 2000 LEVEL = 1, DEPTH 1986 : C 1987 : C ... NUMBER ALL NODES IN THIS LEVEL 1988 : C 1989 : DO 1100 I = 1, LWIDTH 1990 : PPTR = -1 1991 : PTR = WORK (ACTIVE-1) 1992 : IF (NNEXT .EQ. 0) GO TO 1000 1993 : C 1994 : C ... IF NODES REMAIN IN NEXT, FIND THE EARLIEST NODE 1995 : C IN ACTIVE OF MINIMAL DEGREE. 1996 : C 1997 : MINDG = -(N+1) 1998 : DO 400 J = 1, NACTIV 1999 : ASTPTR = ASTART + PTR 2000 : CNODE = WORK (ASTPTR) 2001 : IF ( LVLNUM (CNODE) .EQ. LOWDG ) GO TO 500 2002 : IF ( LVLNUM (CNODE) .LE. MINDG ) GO TO 300 2003 : MPPTR = PPTR 2004 : MPTR = PTR 2005 : MINDG = LVLNUM (CNODE) 2006 : 300 PPTR = PTR 2007 : ACTPTR = ACTIVE + PTR 2008 : PTR = WORK (ACTPTR) 2009 : 400 CONTINUE 2010 : C 2011 : C ... ESTABLISH PTR AS FIRST MIN DEGREE NODE 2012 : C PPTR AS PREDECESSOR IN LIST. 2013 : C 2014 : PTR = MPTR 2015 : PPTR = MPPTR 2016 : C 2017 : 500 ASTPTR = ASTART + PTR 2018 : CNODE = WORK (ASTPTR) 2019 : LOWDG = LVLNUM (CNODE) 2020 : LVLNUM (CNODE) = 0 2021 : JPTR = RSTART (CNODE) 2022 : C 2023 : C ... UPDATE CONNECTION COUNTS FOR ALL NODES WHICH 2024 : C CONNECT TO CNODE'S NEIGHBORS IN NEXT. 2025 : C 2026 : CDGREE = DEGREE(CNODE) 2027 : STRTIC = QUEUEB 2028 : C 2029 : DO 700 J = 1, CDGREE 2030 : JNODE = CONNEC (JPTR) 2031 : JPTR = JPTR + 1 2032 : IF (LVLNUM (JNODE) .NE. LEVEL+1 ) GO TO 700 2033 : IF (QUEUEB .LT. MXQUE) GO TO 5000 2034 : WORK (QUEUEB) = JNODE 2035 : QUEUEB = QUEUEB - 1 2036 : NNEXT = NNEXT - 1 2037 : LVLNUM (JNODE) = 0 2038 : IF (NACTIV .EQ. 1) GO TO 700 2039 : KPTR = RSTART (JNODE) 2040 : JDGREE = DEGREE (JNODE) 2041 : DO 600 K = 1, JDGREE 2042 : KNODE = CONNEC (KPTR) 2043 : KPTR = KPTR + 1 2044 : IF (LVLNUM (KNODE) .GE. 0) GO TO 600 2045 : LVLNUM (KNODE) = LVLNUM (KNODE) + 1 2046 : IF (LOWDG .LT. LVLNUM(KNODE)) 2047 : 1 LOWDG = LVLNUM(KNODE) 2048 : 600 CONTINUE 2049 : 700 CONTINUE 2050 : C 2051 : C ... TO MIMIC THE ALGORITHM AS IMPLEMENTED BY GIBBS, 2052 : C SORT THE NODES JUST ADDED TO THE QUEUE INTO 2053 : C INCREASING ORDER OF ORIGINAL INDEX. (BUT, BECAUSE 2054 : C THE QUEUE IS STORED BACKWARDS IN MEMORY, THE SORT 2054 :
Warning: preg_replace(): Compilation failed: unmatched parentheses at offset 21 in C:\xampp\htdocs\cast3m_kit_local_fr_en_2014\sources.php on line 147
2055 : C 2056 : C TREAT 0, 1 OR 2 NODES ADDED AS SPECIAL CASES 2057 : C 2058 : ADDED = STRTIC - QUEUEB 2059 : IF (ADDED - 2) 1000, 800, 900 2060 : C 2061 : 800 IF (WORK(STRTIC-1) .GT. WORK(STRTIC)) GO TO 1000 2062 : JNODE = WORK(STRTIC) 2063 : WORK(STRTIC) = WORK(STRTIC-1) 2064 : WORK(STRTIC-1) = JNODE 2065 : GO TO 1000 2066 : C 2066 : 900 CALL GPSKCO (ADDED, WORK(QUEUEB+1), ERROR) 2067 : IF (ERROR .NE. 0) GO TO 5500 2068 : C 2069 : C 2070 : C ... NUMBER THIS NODE AND DELETE IT FROM 'ACTIVE'. 2071 : C MARK IT UNAVAILABLE BY CHANGING SIGN OF DEGREE 2072 : C 2073 : 1000 NACTIV = NACTIV - 1 2074 : ASTPTR = ASTART + PTR 2075 : CNODE = WORK (ASTPTR) 2076 : WORK (NXTNUM) = CNODE 2077 : DEGREE (CNODE) = -DEGREE (CNODE) 2078 : NXTNUM = NXTNUM + 1 2079 : C 2080 : C ... DELETE LINK TO THIS NODE FROM LIST 2081 : C 2082 : ACPPTR = ACTIVE + PPTR 2083 : ACTPTR = ACTIVE + PTR 2084 : WORK (ACPPTR) = WORK (ACTPTR) 2085 : 1100 CONTINUE 2086 : C 2087 : C ... NOW MOVE THE QUEUE 'QUEUED' FORWARD, AT THE SAME 2088 : C TIME COMPUTING CONNECTION COUNTS FOR ITS ELEMENTS. 2089 : C THEN DO THE SAME FOR THE REMAINING NODES IN 'NEXT'. 2090 : C 2091 : UNUSED = MIN0 (UNUSED, QUEUEB - MXQUE) 2092 : IF ( NXTNUM .NE. ACTIVE-1 ) GO TO 5100 2093 : IF ( LEVEL .EQ. DEPTH ) GO TO 2000 2094 : LSTART = LVLPTR (LEVEL+1) 2095 : LWIDTH = LVLPTR (LEVEL+2) - LSTART 2096 : ACTIVE = NXTNUM + LWIDTH + 1 2097 : ASTART = ACTIVE + LWIDTH 2098 : NACTIV = LWIDTH 2099 : MXQUE = ASTART + LWIDTH 2100 : IF ( MXQUE .GT. QUEUEB + 1 ) GO TO 5000 2101 : UNUSED = MIN0 (UNUSED, QUEUEB - MXQUE + 1) 2102 : C 2103 : QCOUNT = QUEUED - QUEUEB 2104 : LOWDG = -N-1 2105 : WORK (ACTIVE-1) = 0 2106 : C 2107 : PTR = LSTART 2108 : DO 1600 I = 1, LWIDTH 2109 : C 2110 : C ... CHOOSE NEXT NODE FROM EITHER 'QUEUED' OR 'NEXT' 2111 : C 2112 : IF (I .GT. QCOUNT ) GO TO 1200 2113 : QUEUEI = QUEUED + 1 - I 2114 : CNODE = WORK (QUEUEI) 2115 : GO TO 1300 2116 : C 2117 : 1200 CNODE = LVLLST (PTR) 2118 : PTR = PTR + 1 2119 : IF ( PTR .GT. LVLPTR(LEVEL+2) ) GO TO 5200 2120 : IF (LVLNUM (CNODE) .GT. 0) GO TO 1300 2121 : GO TO 1200 2122 : C 2123 : 1300 IF ( LEVEL+1 .EQ. DEPTH ) GO TO 1500 2124 : C 2125 : RPTR = RSTART (CNODE) 2126 : NCONNC = 0 2127 : JDGREE = DEGREE (CNODE) 2128 : DO 1400 J = 1, JDGREE 2129 : CONNER = CONNEC (RPTR) 2130 : IF ( LVLNUM (CONNER) .EQ. LEVEL+2 ) 2131 : 1 NCONNC = NCONNC - 1 2132 : RPTR = RPTR + 1 2133 : 1400 CONTINUE 2134 : LVLNUM (CNODE) = NCONNC 2135 : LOWDG = MAX0 (LOWDG, NCONNC) 2136 : C 2137 : C ... ADD CNODE TO NEW 'ACTIVE' QUEUE 2138 : C 2139 : 1500 ACTIVI = ACTIVE + (I - 1) 2140 : ASTRTI = ASTART + (I - 1) 2141 : WORK (ACTIVI) = I 2142 : WORK (ASTRTI) = CNODE 2143 : 1600 CONTINUE 2144 : C 2145 : IF (DEPTH .EQ. LEVEL+1 ) GO TO 1700 2146 : NNEXT = LVLPTR (LEVEL+3) - LVLPTR (LEVEL+2) 2147 : QUEUED = LSTART - 1 + LWIDTH + WRKLEN 2148 : QUEUEB = QUEUED 2149 : GO TO 2000 2150 : C 2151 : 1700 NNEXT = 0 2152 : C 2153 : 2000 CONTINUE 2154 : C 2155 : IF (NXTNUM .NE. NFINAL) GO TO 5300 2156 : SPACE = MAX0 (SPACE, TWRKLN - UNUSED) 2157 : RETURN 2158 : C 2159 : C 2160 : C ------------------------------------------------------------------ 2161 : C 2162 : 5000 SPACE = NACTIV + NNEXT 2163 : ERROR = 160 2164 : RETURN 2165 : C 2166 : 5100 ERROR = 61 2167 : GO TO 5400 2168 : C 2169 : 5200 ERROR = 62 2170 : GO TO 5400 2171 : C 2172 : 5300 ERROR = 63 2173 : C 2174 : 5400 RETURN 2175 : C 2176 : 5500 ERROR = 64 2177 : GO TO 5400 2178 : C 2179 : END 2180 : SUBROUTINE GPSKCL (N, DEGREE, RSTART, CONNEC, INVNUM, NEWNUM, 2181 : 1 OLDNUM, BANDWD, PROFIL, ERROR, SPACE) 2182 : C 2183 : C 2184 : IMPLICIT INTEGER(I-N) 2185 : INTEGER N, RSTART(N), BANDWD, PROFIL, ERROR, SPACE 2186 : C 2187 : CIBM INTEGER *2 DEGREE(N), CONNEC(1), INVNUM(N), NEWNUM(N), OLDNUM(N) 2188 : INTEGER DEGREE(N), CONNEC(1), INVNUM(N), NEWNUM(N), OLDNUM(N) 2189 : C 2190 : C ================================================================== 2191 : C 2192 : C 2193 : C COMPUTE THE BANDWIDTH AND PROFILE FOR THE RENUMBERING GIVEN 2194 : C BY 'INVNUM' AND ALSO FOR THE RENUMBERING GIVEN BY 'OLDNUM'. 2195 : C 'NEWNUM' WILL BE A PERMUTATION VECTOR COPY OF THE NODE 2196 : C LIST 'INVNUM'. 2197 : C 2198 : C ================================================================== 2199 : C 2200 : INTEGER I, J, JPTR, IDGREE, OLDBND, OLDPRO, NEWBND, NEWPRO, 2201 : 1 OLDRWD, NEWRWD, OLDORG, NEWORG, JNODE, INVNMI 2202 : C 2203 : C ------------------------------------------------------------------ 2204 : C 2205 : C ... CREATE NEWNUM AS A PERMUTATION VECTOR 2206 : C 2207 : DO 100 I = 1, N 2208 : INVNMI = INVNUM (I) 2209 : NEWNUM (INVNMI) = I 2210 : 100 CONTINUE 2211 : C 2212 : C ... COMPUTE PROFILE AND BANDWIDTH FOR BOTH THE OLD AND THE NEW 2213 : C ORDERINGS. 2214 : C 2215 : OLDBND = 0 2216 : OLDPRO = 0 2217 : NEWBND = 0 2218 : NEWPRO = 0 2219 : C 2220 : DO 300 I = 1, N 2221 : IF (DEGREE(I) .EQ. 0) GO TO 300 2222 : IF (DEGREE(I) .GT. 0) GO TO 6000 2223 : IDGREE = -DEGREE(I) 2224 : DEGREE(I) = IDGREE 2225 : NEWORG = NEWNUM(I) 2226 : OLDORG = OLDNUM(I) 2227 : NEWRWD = 0 2228 : OLDRWD = 0 2229 : JPTR = RSTART (I) 2230 : C 2231 : C ... FIND NEIGHBOR WHICH IS NUMBERED FARTHEST AHEAD OF THE 2232 : C CURRENT NODE. 2233 : C 2234 : DO 200 J = 1, IDGREE 2235 : JNODE = CONNEC(JPTR) 2236 : JPTR = JPTR + 1 2237 : NEWRWD = MAX0 (NEWRWD, NEWORG - NEWNUM(JNODE)) 2238 : OLDRWD = MAX0 (OLDRWD, OLDORG - OLDNUM(JNODE)) 2239 : 200 CONTINUE 2240 : C 2241 : NEWPRO = NEWPRO + NEWRWD 2242 : NEWBND = MAX0 (NEWBND, NEWRWD) 2243 : OLDPRO = OLDPRO + OLDRWD 2244 : OLDBND = MAX0 (OLDBND, OLDRWD) 2245 : 300 CONTINUE 2246 : C 2247 : C ... IF NEW ORDERING HAS BETTER BANDWIDTH THAN OLD ORDERING, 2248 : C REPLACE OLD ORDERING BY NEW ORDERING 2249 : C 2250 : IF (NEWBND .GT. OLDBND) GO TO 500 2251 : BANDWD = NEWBND 2252 : PROFIL = NEWPRO 2253 : DO 400 I = 1, N 2254 : OLDNUM(I) = NEWNUM(I) 2255 : 400 CONTINUE 2256 : GO TO 600 2257 : C 2258 : C ... RETAIN OLD ORDERING 2259 : C 2260 : 500 BANDWD = OLDBND 2261 : PROFIL = OLDPRO 2262 : C 2263 : 600 RETURN 2264 : C 2265 : C ------------------------------------------------------------------ 2266 : C 2267 : 6000 SPACE = -1 2268 : ERROR = 70 2269 : RETURN 2270 : C 2271 : END 2272 : SUBROUTINE GPSKCM (N, DEGREE, RSTART, CONNEC, INVNUM, NEWNUM, 2273 : 1 OLDNUM, BANDWD, PROFIL, ERROR, SPACE) 2274 : C 2275 : C 2276 : IMPLICIT INTEGER(I-N) 2277 : INTEGER N, RSTART(N), BANDWD, PROFIL, ERROR, SPACE 2278 : C 2279 : CIBM INTEGER *2 DEGREE(N), CONNEC(N), INVNUM(N), NEWNUM(N), OLDNUM(N) 2280 : INTEGER DEGREE(N), CONNEC(N), INVNUM(N), NEWNUM(N), OLDNUM(N) 2281 : C 2282 : C ================================================================== 2283 : C 2284 : C 2285 : C COMPUTE THE BANDWIDTH AND PROFILE FOR THE RENUMBERING GIVEN 2286 : C BY 'INVNUM', BY THE REVERSE OF NUMBERING 'INVNUM', AND ALSO 2287 : C BY THE RENUMBERING GIVEN IN 'OLDNUM'. 2288 : C 'NEWNUM' WILL BE A PERMUTATION VECTOR COPY OF THE NODE 2289 : C LIST 'INVNUM'. 2290 : C 2291 : C ================================================================== 2292 : C 2293 : INTEGER I, J, JPTR, IDGREE, OLDBND, OLDPRO, NEWBND, NEWPRO, 2294 : 1 OLDRWD, NEWRWD, OLDORG, NEWORG, JNODE, NRVBND, NRVPRO, 2295 : 2 NRVORG, NRVRWD, INVNMI, NMIP1 2296 : C 2297 : C ------------------------------------------------------------------ 2298 : C 2299 : C ... CREATE NEWNUM AS A PERMUTATION VECTOR 2300 : C 2301 : DO 100 I = 1, N 2302 : INVNMI = INVNUM (I) 2303 : NEWNUM (INVNMI) = I 2304 : 100 CONTINUE 2305 : C 2306 : C ... COMPUTE PROFILE AND BANDWIDTH FOR BOTH THE OLD AND THE NEW 2307 : C ORDERINGS. 2308 : C 2309 : OLDBND = 0 2310 : OLDPRO = 0 2311 : NEWBND = 0 2312 : NEWPRO = 0 2313 : NRVBND = 0 2314 : NRVPRO = 0 2315 : C 2316 : DO 300 I = 1, N 2317 : IF (DEGREE(I) .EQ. 0) GO TO 300 2318 : IF (DEGREE(I) .GT. 0) GO TO 6000 2319 : IDGREE = -DEGREE(I) 2320 : DEGREE(I) = IDGREE 2321 : NEWRWD = 0 2322 : OLDRWD = 0 2323 : NRVRWD = 0 2324 : NEWORG = NEWNUM(I) 2325 : OLDORG = OLDNUM(I) 2326 : NRVORG = N - NEWNUM(I) + 1 2327 : JPTR = RSTART (I) 2328 : C 2329 : C ... FIND NEIGHBOR WHICH IS NUMBERED FARTHEST AHEAD OF THE 2330 : C CURRENT NODE. 2331 : C 2332 : DO 200 J = 1, IDGREE 2333 : JNODE = CONNEC(JPTR) 2334 : JPTR = JPTR + 1 2335 : NEWRWD = MAX0 (NEWRWD, NEWORG - NEWNUM(JNODE)) 2336 : OLDRWD = MAX0 (OLDRWD, OLDORG - OLDNUM(JNODE)) 2337 : NRVRWD = MAX0 (NRVRWD, NRVORG - N + NEWNUM(JNODE) - 1) 2338 : 200 CONTINUE 2339 : C 2340 : NEWPRO = NEWPRO + NEWRWD 2341 : NEWBND = MAX0 (NEWBND, NEWRWD) 2342 : NRVPRO = NRVPRO + NRVRWD 2343 : NRVBND = MAX0 (NRVBND, NRVRWD) 2344 : OLDPRO = OLDPRO + OLDRWD 2345 : OLDBND = MAX0 (OLDBND, OLDRWD) 2346 : 300 CONTINUE 2347 : C 2348 : C ... IF NEW ORDERING HAS BETTER BANDWIDTH THAN OLD ORDERING, 2349 : C REPLACE OLD ORDERING BY NEW ORDERING 2350 : C 2351 : IF ((NEWPRO .GT. OLDPRO) .OR. (NEWPRO .GT. NRVPRO)) GO TO 500 2352 : BANDWD = NEWBND 2353 : PROFIL = NEWPRO 2354 : DO 400 I = 1, N 2355 : OLDNUM(I) = NEWNUM(I) 2356 : 400 CONTINUE 2357 : GO TO 800 2358 : C 2359 : C ... CHECK NEW REVERSED ORDERING FOR BEST PROFILE 2360 : C 2361 : 500 IF (NRVPRO .GT. OLDPRO) GO TO 700 2362 : BANDWD = NRVBND 2363 : PROFIL = NRVPRO 2364 : DO 600 I = 1, N 2365 : OLDNUM(I) = N - NEWNUM(I) + 1 2366 : IF (I .GT. N/2) GO TO 600 2367 : J = INVNUM(I) 2368 : NMIP1 = (N + 1) - I 2369 : INVNUM(I) = INVNUM (NMIP1) 2370 : INVNUM (NMIP1) = J 2371 : 600 CONTINUE 2372 : GO TO 800 2373 : C 2374 : C 2375 : C ... RETAIN OLD ORDERING 2376 : C 2377 : 700 BANDWD = OLDBND 2378 : PROFIL = OLDPRO 2379 : C 2380 : 800 RETURN 2381 : C 2382 : C ------------------------------------------------------------------ 2383 : C 2384 : 6000 ERROR = 71 2385 : SPACE = -1 2386 : RETURN 2387 : C 2388 : END 2389 : SUBROUTINE GPSKCN (N, KEY, DATA, ERROR) 2390 : C 2391 : C ================================================================== 2392 : C 2393 : C I N S E R T I O N S O R T 2394 : C 2395 : C INPUT: 2396 : C N -- NUMBER OF ELEMENTS TO BE SORTED 2397 : C KEY -- AN ARRAY OF LENGTH N CONTAINING THE VALUES 2398 : C WHICH ARE TO BE SORTED 2399 : C DATA -- A SECOND ARRAY OF LENGTH N CONTAINING DATA 2400 : C ASSOCIATED WITH THE INDIVIDUAL KEYS. 2401 : C 2402 : C OUTPUT: 2403 : C KEY -- WILL BE ARRANGED SO THAT VALUES ARE IN DECREASING 2404 : C ORDER 2405 : C DATA -- REARRANGED TO CORRESPOND TO REARRANGED KEYS 2406 : C ERROR -- WILL BE ZERO UNLESS THE PROGRAM IS MALFUNCTIONING, 2407 : C IN WHICH CASE IT WILL BE EQUAL TO 1. 2408 : C 2409 : C 2410 : C ================================================================== 2411 : C 2412 : IMPLICIT INTEGER(I-N) 2413 : INTEGER N, ERROR 2414 : C 2415 : CIBM INTEGER *2 KEY(N), DATA(N) 2416 : INTEGER KEY(N), DATA(N) 2417 : C 2418 : C ------------------------------------------------------------------ 2419 : C 2420 : INTEGER I, J, D, K, IP1, JM1 2421 : C 2422 : C ------------------------------------------------------------------ 2423 : C 2424 : IF (N .EQ. 1) RETURN 2425 : IF (N .LE. 0) GO TO 6000 2426 : C 2427 : ERROR = 0 2428 : C 2429 : C ... INSERTION SORT ... FOR I := N-1 STEP -1 TO 1 DO ... 2430 : C 2431 : 2400 I = N - 1 2432 : IP1 = N 2433 : C 2434 : 2500 IF ( KEY (I) .GE. KEY (IP1) ) GO TO 2800 2435 : C 2436 : C ... OUT OF ORDER ... MOVE UP TO CORRECT PLACE 2437 : C 2438 : K = KEY (I) 2439 : D = DATA (I) 2440 : J = IP1 2441 : JM1 = I 2442 : C 2443 : C ... REPEAT ... UNTIL 'CORRECT PLACE FOR K FOUND' 2444 : C 2445 : 2600 KEY (JM1) = KEY (J) 2446 : DATA (JM1) = DATA (J) 2447 : JM1 = J 2448 : J = J + 1 2449 : IF (J .GT. N) GO TO 2700 2450 : IF (KEY (J) .GT. K) GO TO 2600 2451 : C 2452 : 2700 KEY (JM1) = K 2453 : DATA (JM1) = D 2454 : C 2455 : 2800 IP1 = I 2456 : I = I - 1 2457 : IF ( I .GT. 0 ) GO TO 2500 2458 : C 2459 : 3000 RETURN 2460 : C 2461 : 6000 ERROR = 1 2462 : GO TO 3000 2463 : C 2464 : END 2465 : SUBROUTINE GPSKCO (N, KEY, ERROR) 2466 : C 2467 : C ================================================================== 2468 : C 2469 : C I N S E R T I O N S O R T 2470 : C 2471 : C INPUT: 2472 : C N -- NUMBER OF ELEMENTS TO BE SORTED 2473 : C KEY -- AN ARRAY OF LENGTH N CONTAINING THE VALUES 2474 : C WHICH ARE TO BE SORTED 2475 : C 2476 : C OUTPUT: 2477 : C KEY -- WILL BE ARRANGED SO THAT VALUES ARE IN DECREASING 2478 : C ORDER 2479 : C 2480 : C ================================================================== 2481 : C 2482 : IMPLICIT INTEGER(I-N) 2483 : INTEGER N, ERROR 2484 : C 2485 : CIBM INTEGER *2 KEY(N) 2486 : INTEGER KEY(N) 2487 : C 2488 : C ------------------------------------------------------------------ 2489 : C 2490 : INTEGER I, J, K, IP1, JM1 2491 : C 2492 : C ------------------------------------------------------------------ 2493 : C 2494 : IF (N .EQ. 1) RETURN 2495 : IF (N .LE. 0) GO TO 6000 2496 : C 2497 : ERROR = 0 2498 : C 2499 : C ... INSERTION SORT ... FOR I := N-1 STEP -1 TO 1 DO ... 2500 : C 2501 : 2400 I = N - 1 2502 : IP1 = N 2503 : C 2504 : 2500 IF ( KEY (I) .GE. KEY (IP1) ) GO TO 2800 2505 : C 2506 : C ... OUT OF ORDER ... MOVE UP TO CORRECT PLACE 2507 : C 2508 : K = KEY (I) 2509 : J = IP1 2510 : JM1 = I 2511 : C 2512 : C ... REPEAT ... UNTIL 'CORRECT PLACE FOR K FOUND' 2513 : C 2514 : 2600 KEY (JM1) = KEY (J) 2515 : JM1 = J 2516 : J = J + 1 2517 : IF (J .GT. N) GO TO 2700 2518 : IF (KEY (J) .GT. K) GO TO 2600 2519 : C 2520 : 2700 KEY (JM1) = K 2521 : C 2522 : 2800 IP1 = I 2523 : I = I - 1 2524 : IF ( I .GT. 0 ) GO TO 2500 2525 : C 2526 : 3000 RETURN 2527 : C 2528 : 6000 ERROR = 1 2529 : GO TO 3000 2530 : C 2531 : END 2532 : SUBROUTINE GPSKCP (N, INDEX, NVEC, DEGREE, ERROR) 2533 : C 2534 : C ================================================================== 2535 : C 2536 : C I N S E R T I O N S O R T 2537 : C 2538 : C INPUT: 2539 : C N -- NUMBER OF ELEMENTS TO BE SORTED 2540 : C INDEX -- AN ARRAY OF LENGTH N CONTAINING THE INDICES 2541 : C WHOSE DEGREES ARE TO BE SORTED 2542 : C DEGREE -- AN NVEC VECTOR, GIVING THE DEGREES OF NODES 2543 : C WHICH ARE TO BE SORTED. 2544 : C 2545 : C OUTPUT: 2546 : C INDEX -- WILL BE ARRANGED SO THAT VALUES ARE IN DECREASING 2547 : C ORDER 2548 : C ERROR -- WILL BE ZERO UNLESS THE PROGRAM IS MALFUNCTIONING, 2549 : C IN WHICH CASE IT WILL BE EQUAL TO 1. 2550 : C 2551 : C ================================================================== 2552 : C 2553 : IMPLICIT INTEGER(I-N) 2554 : INTEGER N, NVEC, ERROR 2555 : C 2556 : CIBM INTEGER *2 INDEX(N), DEGREE(NVEC) 2557 : INTEGER INDEX(N), DEGREE(NVEC) 2558 : C 2559 : C ------------------------------------------------------------------ 2560 : C 2561 : INTEGER I, J, V, IP1, JM1, INDEXI, INDXI1, INDEXJ 2562 : C 2563 : C ------------------------------------------------------------------ 2564 : C 2565 : IF (N .EQ. 1) RETURN 2566 : IF (N .LE. 0) GO TO 6000 2567 : C 2568 : ERROR = 0 2569 : C 2570 : C ------------------------------------------------------------------ 2571 : C INSERTION SORT THE ENTIRE FILE 2572 : C ------------------------------------------------------------------ 2573 : C 2574 : C 2575 : C ... INSERTION SORT ... FOR I := N-1 STEP -1 TO 1 DO ... 2576 : C 2577 : 2400 I = N - 1 2578 : IP1 = N 2579 : C 2580 : 2500 INDEXI = INDEX (I) 2581 : INDXI1 = INDEX (IP1) 2582 : IF ( DEGREE(INDEXI) .GE. DEGREE(INDXI1) ) GO TO 2800 2583 : C 2584 : C ... OUT OF ORDER ... MOVE UP TO CORRECT PLACE 2585 : C 2586 : V = DEGREE (INDEXI) 2587 : J = IP1 2588 : JM1 = I 2589 : INDEXJ = INDEX (J) 2590 : C 2591 : C ... REPEAT ... UNTIL 'CORRECT PLACE FOR V FOUND' 2592 : C 2593 : 2600 INDEX (JM1) = INDEXJ 2594 : JM1 = J 2595 : J = J + 1 2596 : IF (J .GT. N) GO TO 2700 2597 : INDEXJ = INDEX (J) 2598 : IF (DEGREE(INDEXJ) .GT. V) GO TO 2600 2599 : C 2600 : 2700 INDEX (JM1) = INDEXI 2601 : C 2602 : 2800 IP1 = I 2603 : I = I - 1 2604 : IF ( I .GT. 0 ) GO TO 2500 2605 : C 2606 : 3000 RETURN 2607 : C 2608 : 6000 ERROR = 1 2609 : GO TO 3000 2610 : C 2611 : END 2612 : SUBROUTINE GPSKCQ (N, INDEX, NVEC, DEGREE, ERROR) 2613 : C 2614 : C ================================================================== 2615 : C 2616 : C I N S E R T I O N S O R T 2617 : C 2618 : C INPUT: 2619 : C N -- NUMBER OF ELEMENTS TO BE SORTED 2620 : C INDEX -- AN ARRAY OF LENGTH N CONTAINING THE INDICES 2621 : C WHOSE DEGREES ARE TO BE SORTED 2622 : C DEGREE -- AN NVEC VECTOR, GIVING THE DEGREES OF NODES 2623 : C WHICH ARE TO BE SORTED. 2624 : C 2625 : C OUTPUT: 2626 : C INDEX -- WILL BE ARRANGED SO THAT VALUES ARE IN INCREASING 2627 : C ORDER 2628 : C ERROR -- WILL BE ZERO UNLESS THE PROGRAM IS MALFUNCTIONING, 2629 : C IN WHICH CASE IT WILL BE EQUAL TO 1. 2630 : C 2631 : C ================================================================== 2632 : C 2633 : IMPLICIT INTEGER(I-N) 2634 : INTEGER N, NVEC, ERROR 2635 : C 2636 : CIBM INTEGER *2 INDEX(N), DEGREE(NVEC) 2637 : INTEGER INDEX(N), DEGREE(NVEC) 2638 : C 2639 : C ------------------------------------------------------------------ 2640 : C 2641 : INTEGER I, J, V, INDEXI, INDXI1, INDEXJ, IP1, JM1 2642 : C 2643 : C ------------------------------------------------------------------ 2644 : C 2645 : IF (N .EQ. 1) RETURN 2646 : IF (N .LE. 0) GO TO 6000 2647 : C 2648 : ERROR = 0 2649 : C 2650 : C ------------------------------------------------------------------ 2651 : C INSERTION SORT THE ENTIRE FILE 2652 : C ------------------------------------------------------------------ 2653 : C 2654 : C 2655 : C ... INSERTION SORT ... FOR I := N-1 STEP -1 TO 1 DO ... 2656 : C 2657 : 2400 I = N - 1 2658 : IP1 = N 2659 : C 2660 : 2500 INDEXI = INDEX (I) 2661 : INDXI1 = INDEX (IP1) 2662 : IF ( DEGREE(INDEXI) .LE. DEGREE(INDXI1) ) GO TO 2800 2663 : C 2664 : C ... OUT OF ORDER ... MOVE UP TO CORRECT PLACE 2665 : C 2666 : V = DEGREE (INDEXI) 2667 : J = IP1 2668 : JM1 = I 2669 : INDEXJ = INDEX (J) 2670 : C 2671 : C ... REPEAT ... UNTIL 'CORRECT PLACE FOR V FOUND' 2672 : C 2673 : 2600 INDEX (JM1) = INDEXJ 2674 : JM1 = J 2675 : J = J + 1 2676 : IF (J .GT. N) GO TO 2700 2677 : INDEXJ = INDEX (J) 2678 : IF (DEGREE(INDEXJ) .LT. V) GO TO 2600 2679 : C 2680 : 2700 INDEX (JM1) = INDEXI 2681 : C 2682 : 2800 IP1 = I 2683 : I = I - 1 2684 : IF ( I .GT. 0 ) GO TO 2500 2685 : C 2686 : 3000 RETURN 2687 : C 2688 : 6000 ERROR = 1 2689 : GO TO 3000 2690 : C 2691 : END 2692 : 2693 : 2694 :
© Cast3M 2003 - All rights reserved.
Disclaimer