Télécharger gpskc.eso

Retour à la liste

Numérotation des lignes :

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

© Cast3M 2003 - Tous droits réservés.
Mentions légales