Context Navigation

source: trunk/SOURCES/reduc.f @ 23

Last change on this file since 23 was 4, checked in by dumas, 10 years ago
initial import GRISLI trunk
File size: 90.1 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: