C ALGORITHM 636 COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.11, NO. 4, C DEC., 1985, P. 378. C ********** C C THIS IS A TEST PROGRAM FOR SUBROUTINES DSSM AND FDHS. C THE TEST DATA REPRESENTS MINIMAL SURFACE PROBLEMS. C C ********** INTEGER MAXNNZ,MAXN,LIWA PARAMETER (MAXNNZ = 8000, MAXN = 1600, LIWA = 9600) C C LOGICAL OUTPUT UNIT IS ASSUMED TO BE NUMBER 6. C INTEGER NWRITE PARAMETER (NWRITE = 6) C INTEGER I,IX,INFO,J,JY,JP,L,MAXGRP,MAXROW,METHOD,MINGRP, * N,NNZ,NUMGRP INTEGER INDCOL(MAXNNZ),INDROW(MAXNNZ),IPNTR(MAXN+1), * JPNTR(MAXN+1),LISTP(MAXN),NGRP(MAXN),IWA(LIWA) REAL ABSERR,DNSM,ENTRY1,ENTRY2,ERRIJ,FHEST,FM,RELERR REAL ETA(MAXN),FHES(MAXNNZ),FHESD(MAXN),GVEC(MAXN), * X(MAXN),XD(MAXN) REAL HS C C TEST FOR DSSM AND FDHS. C WRITE (NWRITE,1000) DO 100 L = 10, 40, 10 N = L**2 DO 90 METHOD = 1, 2 C C DEFINITION OF SPARSITY PATTERN. C NNZ = 0 DO 10 J = 1, N NNZ = NNZ + 1 INDCOL(NNZ) = J INDROW(NNZ) = J IF (MOD(J,L) .NE. 0) THEN NNZ = NNZ + 1 INDCOL(NNZ) = J INDROW(NNZ) = J + 1 END IF IF (J+L .LE. N) THEN NNZ = NNZ + 1 INDCOL(NNZ) = J INDROW(NNZ) = J + L IF (MOD(J,L) .NE. 1) THEN NNZ = NNZ + 1 INDCOL(NNZ) = J INDROW(NNZ) = J + L - 1 END IF IF (MOD(J,L) .NE. 0) THEN NNZ = NNZ + 1 INDCOL(NNZ) = J INDROW(NNZ) = J + L + 1 END IF END IF 10 CONTINUE C C CALL DSSM. C CALL DSSM(N,NNZ,INDROW,INDCOL,METHOD,LISTP,NGRP,MAXGRP, * MINGRP,INFO,IPNTR,JPNTR,IWA,LIWA) IF (INFO .LE. 0) THEN WRITE (NWRITE,2000) INFO STOP END IF C C STATISTICS FOR THE MATRIX. C MAXROW = 0 DO 20 I = 1, N MAXROW = MAX(MAXROW,IPNTR(I+1)-IPNTR(I)) 20 CONTINUE DNSM = FLOAT(100*NNZ)/FLOAT(N*N) WRITE (NWRITE,3000) N,NNZ,DNSM,MAXROW,MINGRP,MAXGRP C C TEST FOR FDHS. C FM = (L+1)**2 DO 30 J = 1, N JY = (J-1)/L + 1 IX = J - L*(JY - 1) X(J) = (JY**2 + IX**2)/FM 30 CONTINUE CALL FCN(L,X,GVEC) C C APPROXIMATE THE HESSIAN MATRIX. C DO 60 NUMGRP = 1, MAXGRP DO 40 J = 1, N ETA(J) = 1.0E-4 XD(J) = X(J) IF (NGRP(J) .EQ. NUMGRP) XD(J) = X(J) + ETA(J) 40 CONTINUE CALL FCN(L,XD,FHESD) DO 50 I = 1, N FHESD(I) = FHESD(I) - GVEC(I) 50 CONTINUE CALL FDHS(N,INDROW,JPNTR,INDCOL,IPNTR, * LISTP,NGRP,MAXGRP,NUMGRP,ETA,FHESD,FHES,IWA) 60 CONTINUE C C TEST THE APPROXIMATION TO THE HESSIAN. C ABSERR = 0.0 RELERR = 0.0 ENTRY1 = 0.0 ENTRY2 = 0.0 DO 80 J = 1, N JY = (J - 1)/L IX = J - L*JY - 1 DO 70 JP = JPNTR(J), JPNTR(J+1)-1 I = INDROW(JP) IF (I .EQ. J) THEN FHEST = HS(IX,JY,L,X,1) + HS(IX+1,JY,L,X,2) + * HS(IX,JY+1,L,X,2) + HS(IX+1,JY+1,L,X,1) ELSE IF (I .EQ. J+1) THEN FHEST = HS(IX+1,JY,L,X,3) + HS(IX+1,JY+1,L,X,3) ELSE IF (I .EQ. J+L) THEN FHEST = - HS(IX,JY+1,L,X,3) - HS(IX+1,JY+1,L,X,3) ELSE IF (I .EQ. J+L-1) THEN FHEST = -HS(IX,JY+1,L,X,2) ELSE IF (I .EQ. J+L+1) THEN FHEST = -HS(IX+1,JY+1,L,X,1) END IF ERRIJ = FHES(JP) - FHEST IF (ABS(ERRIJ) .GT. ABSERR) THEN ENTRY1 = FHEST ABSERR = ABS(ERRIJ) END IF IF (ABS(ERRIJ) .GT. RELERR*ABS(FHEST)) THEN ENTRY2 = FHEST IF (FHEST .NE. 0.0) * RELERR = ABS(ERRIJ)/ABS(FHEST) END IF 70 CONTINUE 80 CONTINUE WRITE (NWRITE,4000) ABSERR,ENTRY1,RELERR,ENTRY2 90 CONTINUE 100 CONTINUE C C FORMAT STATEMENTS. C 1000 FORMAT (// ' TESTS FOR DSSM - MINIMAL SURFACE PROBLEMS' // * ' STATISTICS GENERATED '// * ' N - NUMBER OF COLUMNS '/ * ' NNZ - NUMBER OF NONZERO ELEMENTS'/ * ' DNSM - MATRIX DENSITY (PERCENTAGE)'/ * ' MAXROW - MAXIMUM NUMBER OF NONZEROES IN ANY ROW'/ * ' MINGRP - LOWER BOUND ON NUMBER OF GROUPS'/ * ' MAXGRP - NUMBER OF GROUPS DETERMINED BY DSSM'///) 2000 FORMAT (// ' *** MISTAKE IN INPUT DATA, INFO IS ***',I6) 3000 FORMAT (// 3X,'N',6X,'NNZ',5X,'DNSM',5X, * 'MAXROW',4X,'MINGRP',4X,'MAXGRP'// * 2(I5,3X),F6.2,4X,3(I5,5X)) 4000 FORMAT (// ' LARGEST ABSOLUTE ERROR ',E10.2,5X,'ENTRY',E10.2, * / ' LARGEST RELATIVE ERROR ',E10.2,5X,'ENTRY',E10.2) STOP C C LAST CARD OF MAIN PROGRAM. C END SUBROUTINE FCN(L,X,GVEC) INTEGER L REAL X(L**2),GVEC(L**2) C C GRADIENT SUBROUTINE FOR TESTING FDHS. C INTEGER I,IX,JY,K REAL A,B,C,D,FM,PHI,PHI2 C DO 10 I = 1, L**2 GVEC(I) = 0.0 10 CONTINUE C C THE GRADIENT IS COMPUTED AS THE SUM OF THE C GRADIENTS OF THE (L+1)**2 ELEMENT FUNCTIONS. C DO 30 JY = 0, L DO 20 IX = 0, L K = L*(JY - 1) + IX C C COMPUTE THE GRADIENT OF THE FUNCTION ASSOCIATED WITH THE C ELEMENT WITH CORNERS AT X(K), X(K+1), X(K+L), X(K+L+1). C FM = (L+1)**2 IF (IX .EQ. 0 .OR. JY .EQ. 0) THEN A = (IX**2 + JY**2)/FM ELSE A = X(K) END IF IF (IX .EQ. L .OR. JY .EQ. 0) THEN B = ((IX+1)**2 + JY**2)/FM ELSE B = X(K+1) END IF IF (IX .EQ. 0 .OR. JY .EQ. L) THEN C = (IX**2 + (JY+1)**2)/FM ELSE C = X(K+L) END IF IF (IX .EQ. L .OR. JY .EQ. L) THEN D = ((IX+1)**2 + (JY+1)**2)/FM ELSE D = X(K+L+1) END IF PHI = 1.0 + 0.5*FM*((A-D)**2 + (B-C)**2) PHI2 = 2*SQRT(PHI) IF (IX .NE. 0 .AND. JY .NE. 0) * GVEC(K) = GVEC(K) + (A-D)/PHI2 IF (IX .NE. L .AND. JY .NE. 0) * GVEC(K+1) = GVEC(K+1) + (B-C)/PHI2 IF (IX .NE. 0 .AND. JY .NE. L) * GVEC(K+L) = GVEC(K+L) + (C-B)/PHI2 IF (IX .NE. L .AND. JY .NE. L) * GVEC(K+L+1) = GVEC(K+L+1) + (D-A)/PHI2 20 CONTINUE 30 CONTINUE RETURN C C LAST CARD OF SUBROUTINE FCN. C END REAL FUNCTION HS(IX,JY,L,X,MODE) INTEGER IX,JY,L,MODE REAL X(L**2) C C THIS FUNCTION COMPUTES THE HESSIAN OF THE FUNCTION C ASSOCIATED WITH THE ELEMENT WITH CORNERS AT X(K), C X(K+1), X(K+L), X(K+L+1), WHERE K = L*(JY - 1) + IX. C INTEGER K REAL A,B,C,D,FM,PHI,PHI2 FM = (L+1)**2 K = L*(JY - 1) + IX IF (IX .EQ. 0 .OR. JY .EQ. 0) THEN A = (IX**2 + JY**2)/FM ELSE A = X(K) END IF IF (IX .EQ. L .OR. JY .EQ. 0) THEN B = ((IX+1)**2 + JY**2)/FM ELSE B = X(K+1) END IF IF (IX .EQ. 0 .OR. JY .EQ. L) THEN C = (IX**2 + (JY+1)**2)/FM ELSE C = X(K+L) END IF IF (IX .EQ. L .OR. JY .EQ. L) THEN D = ((IX+1)**2 + (JY+1)**2)/FM ELSE D = X(K+L+1) END IF PHI = 1.0 + 0.5*FM*((A-D)**2 + (B-C)**2) PHI2 = 2*SQRT(PHI) IF (MODE .EQ. 1) HS = ((1.0 + 0.5*FM*(B-C)**2)/PHI)/PHI2 IF (MODE .EQ. 2) HS = ((1.0 + 0.5*FM*(A-D)**2)/PHI)/PHI2 IF (MODE .EQ. 3) HS = -0.5*FM*(((A-D)*(B-C)/PHI)/PHI2) RETURN C C LAST CARD OF FUNCTION HS. C END SUBROUTINE DSSM(N,NPAIRS,INDROW,INDCOL,METHOD,LISTP,NGRP, * MAXGRP,MINGRP,INFO,IPNTR,JPNTR,IWA,LIWA) INTEGER N,NPAIRS,METHOD,MAXGRP,MINGRP,INFO,LIWA INTEGER INDROW(NPAIRS),INDCOL(NPAIRS),LISTP(N),NGRP(N), * IPNTR(N+1),JPNTR(N+1),IWA(LIWA) C ********** C C SUBROUTINE DSSM C C GIVEN THE SPARSITY PATTERN OF A SYMMETRIC MATRIX A OF ORDER N, C THIS SUBROUTINE DETERMINES A SYMMETRIC PERMUTATION OF A AND A C PARTITION OF THE COLUMNS OF A CONSISTENT WITH THE DETERMINATION C OF A BY A LOWER TRIANGULAR SUBSTITUTION METHOD. C C THE SPARSITY PATTERN OF THE MATRIX A IS SPECIFIED BY THE C ARRAYS INDROW AND INDCOL. ON INPUT THE INDICES FOR THE C NON-ZERO ELEMENTS IN THE LOWER TRIANGULAR PART OF A ARE C C (INDROW(K),INDCOL(K)), K = 1,2,...,NPAIRS. C C THE (INDROW(K),INDCOL(K)) PAIRS MAY BE SPECIFIED IN ANY ORDER. C DUPLICATE INPUT PAIRS ARE PERMITTED, BUT THE SUBROUTINE C ELIMINATES THEM. THE SUBROUTINE REQUIRES THAT ALL THE DIAGONAL C ELEMENTS BE PART OF THE SPARSITY PATTERN AND REPLACES ANY PAIR C (INDROW(K),INDCOL(K)) WHERE INDROW(K) IS LESS THAN INDCOL(K) C BY THE PAIR (INDCOL(K),INDROW(K)). C C THE DIRECT METHOD (METHOD = 1) FIRST DETERMINES A PARTITION C OF THE COLUMNS OF A SUCH THAT TWO COLUMNS IN A GROUP HAVE A C NON-ZERO ELEMENT IN ROW K ONLY IF COLUMN K IS IN AN EARLIER C GROUP. USING THIS PARTITION, THE SUBROUTINE THEN COMPUTES A C SYMMETRIC PERMUTATION OF A CONSISTENT WITH THE DETERMINATION C OF A BY A LOWER TRIANGULAR SUBSTITUTION METHOD. C C THE INDIRECT METHOD FIRST COMPUTES A SYMMETRIC PERMUTATION OF A C WHICH MINIMIZES THE MAXIMUM NUMBER OF NON-ZERO ELEMENTS IN ANY C ROW OF L, WHERE L IS THE LOWER TRIANGULAR PART OF THE PERMUTED C MATRIX. THE SUBROUTINE THEN PARTITIONS THE COLUMNS OF L INTO C GROUPS SUCH THAT COLUMNS OF L IN A GROUP DO NOT HAVE A NON-ZERO C IN THE SAME ROW POSITION. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE DSSM(N,NPAIRS,INDROW,INDCOL,METHOD,LISTP,NGRP, C MAXGRP,MINGRP,INFO,IPNTR,JPNTR,IWA,LIWA) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE ORDER OF A. C C NPAIRS IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF (INDROW,INDCOL) PAIRS USED TO DESCRIBE THE SPARSITY C PATTERN OF A. C C INDROW IS AN INTEGER ARRAY OF LENGTH NPAIRS. ON INPUT INDROW C MUST CONTAIN THE ROW INDICES OF THE NON-ZERO ELEMENTS IN C THE LOWER TRIANGULAR PART OF A. ON OUTPUT INDROW IS C PERMUTED SO THAT THE CORRESPONDING COLUMN INDICES ARE IN C NON-DECREASING ORDER. THE COLUMN INDICES CAN BE RECOVERED C FROM THE ARRAY JPNTR. C C INDCOL IS AN INTEGER ARRAY OF LENGTH NPAIRS. ON INPUT INDCOL C MUST CONTAIN THE COLUMN INDICES OF THE NON-ZERO ELEMENTS C IN THE LOWER TRIANGULAR PART OF A. ON OUTPUT INDCOL IS C PERMUTED SO THAT THE CORRESPONDING ROW INDICES ARE IN C NON-DECREASING ORDER. THE ROW INDICES CAN BE RECOVERED C FROM THE ARRAY IPNTR. C C METHOD IS AN INTEGER INPUT VARIABLE. IF METHOD = 1, THE C DIRECT METHOD IS USED TO DETERMINE THE PARTITION AND C SYMMETRIC PERMUTATION. OTHERWISE, THE INDIRECT METHOD IS C USED TO DETERMINE THE SYMMETRIC PERMUTATION AND PARTITION. C C LISTP IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE SYMMETRIC PERMUTATION OF THE MATRIX A. ELEMENT (I,J) C OF A IS THE (LISTP(I),LISTP(J)) ELEMENT OF THE PERMUTED C MATRIX. C C NGRP IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE PARTITION OF THE COLUMNS OF A. COLUMN J BELONGS TO C GROUP NGRP(J). C C MAXGRP IS AN INTEGER OUTPUT VARIABLE WHICH SPECIFIES THE C NUMBER OF GROUPS IN THE PARTITION OF THE COLUMNS OF A. C C MINGRP IS AN INTEGER OUTPUT VARIABLE WHICH SPECIFIES A LOWER C BOUND FOR THE NUMBER OF GROUPS IN ANY PARTITION OF THE C COLUMNS OF A CONSISTENT WITH THE DETERMINATION OF A BY A C LOWER TRIANGULAR SUBSTITUTION METHOD. C C INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS. FOR C NORMAL TERMINATION INFO = 1. IF N OR NPAIRS IS NOT C POSITIVE OR LIWA IS LESS THAN 6*N, THEN INFO = 0. IF THE C K-TH ELEMENT OF INDROW OR THE K-TH ELEMENT OF INDCOL IS C NOT AN INTEGER BETWEEN 1 AND N, OR IF THE K-TH DIAGONAL C ELEMENT IS NOT IN THE SPARSITY PATTERN, THEN INFO = -K. C C IPNTR IS AN INTEGER OUTPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL. C THE COLUMN INDICES FOR ROW I ARE C C INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1. C C NOTE THAT IPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO C ELEMENTS IN THE LOWER TRIANGULAR PART OF THE MATRIX A. C C JPNTR IS AN INTEGER OUTPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW. C THE ROW INDICES FOR COLUMN J ARE C C INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1. C C NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO C ELEMENTS IN THE LOWER TRIANGULAR PART OF THE MATRIX A. C C IWA IS AN INTEGER WORK ARRAY OF LENGTH LIWA. C C LIWA IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN 6*N. C C SUBPROGRAMS CALLED C C MINPACK-SUPPLIED ... DEGR,IDO,IDOG,NUMSRT,SDPT,SEQ,SETR, C SLO,SLOG,SRTDAT C C FORTRAN-SUPPLIED ... MAX,MIN C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. DECEMBER 1984. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER I,IR,J,JP,K,MAXID,MAXVD,MAXCLQ,NNZ,NUMGRP C C CHECK THE INPUT DATA. C INFO = 0 IF (N .LT. 1 .OR. NPAIRS .LT. 1 .OR. LIWA .LT. 6*N) RETURN DO 10 K = 1, N IWA(K) = 0 10 CONTINUE DO 20 K = 1, NPAIRS INFO = -K IF (INDROW(K) .LT. 1 .OR. INDROW(K) .GT. N .OR. * INDCOL(K) .LT. 1 .OR. INDCOL(K) .GT. N) RETURN IF (INDROW(K) .EQ. INDCOL(K)) IWA(INDROW(K)) = 1 20 CONTINUE DO 30 K = 1, N INFO = -K IF (IWA(K) .NE. 1) RETURN 30 CONTINUE INFO = 1 C C GENERATE THE SPARSITY PATTERN FOR THE LOWER C TRIANGULAR PART OF A. C DO 40 K = 1, NPAIRS I = INDROW(K) J = INDCOL(K) INDROW(K) = MAX(I,J) INDCOL(K) = MIN(I,J) 40 CONTINUE C C SORT THE DATA STRUCTURE BY COLUMNS. C CALL SRTDAT(N,NPAIRS,INDROW,INDCOL,JPNTR,IWA) C C COMPRESS THE DATA AND DETERMINE THE NUMBER OF NON-ZERO C ELEMENTS IN THE LOWER TRIANGULAR PART OF A. C DO 50 I = 1, N IWA(I) = 0 50 CONTINUE NNZ = 0 DO 70 J = 1, N K = NNZ DO 60 JP = JPNTR(J), JPNTR(J+1)-1 IR = INDROW(JP) IF (IWA(IR) .NE. J) THEN NNZ = NNZ + 1 INDROW(NNZ) = IR IWA(IR) = J END IF 60 CONTINUE JPNTR(J) = K + 1 70 CONTINUE JPNTR(N+1) = NNZ + 1 C C EXTEND THE DATA STRUCTURE TO ROWS. C CALL SETR(N,N,INDROW,JPNTR,INDCOL,IPNTR,IWA) C C DETERMINE THE SMALLEST-LAST ORDERING OF THE VERTICES OF THE C ADJACENCY GRAPH OF A, AND FROM IT DETERMINE A LOWER BOUND C FOR THE NUMBER OF GROUPS. C CALL SLOG(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(1),MAXCLQ, * MAXVD,IWA(N+1),IWA(2*N+1),IWA(3*N+1)) MINGRP = 1 + MAXVD C C USE THE SELECTED METHOD. C IF (METHOD .EQ. 1) THEN C C DIRECT METHOD. DETERMINE A PARTITION OF THE COLUMNS C OF A BY THE POWELL-TOINT METHOD. C CALL SDPT(N,INDROW,JPNTR,INDCOL,IPNTR,NGRP,MAXGRP, * IWA(N+1),IWA(2*N+1)) C C DEFINE A SYMMETRIC PERMUTATION OF A ACCORDING TO THE C ORDERING OF THE COLUMN GROUP NUMBERS IN THE PARTITION. C CALL NUMSRT(N,MAXGRP,NGRP,1,IWA(1),IWA(2*N+1),IWA(N+1)) DO 80 I = 1, N LISTP(IWA(I)) = I 80 CONTINUE ELSE C C INDIRECT METHOD. DETERMINE THE INCIDENCE DEGREE ORDERING C OF THE VERTICES OF THE ADJACENCY GRAPH OF A AND, TOGETHER C WITH THE SMALLEST-LAST ORDERING, DEFINE A SYMMETRIC C PERMUTATION OF A. C CALL IDOG(N,INDROW,JPNTR,INDCOL,IPNTR,LISTP,MAXCLQ, * MAXID,IWA(N+1),IWA(2*N+1),IWA(3*N+1)) IF (MAXID .GT. MAXVD) THEN DO 90 I = 1, N LISTP(I) = IWA(I) 90 CONTINUE END IF C C GENERATE THE SPARSITY PATTERN FOR THE LOWER C TRIANGULAR PART L OF THE PERMUTED MATRIX. C DO 110 J = 1, N DO 100 JP = JPNTR(J), JPNTR(J+1)-1 I = INDROW(JP) INDROW(JP) = MAX(LISTP(I),LISTP(J)) INDCOL(JP) = MIN(LISTP(I),LISTP(J)) 100 CONTINUE 110 CONTINUE C C SORT THE DATA STRUCTURE BY COLUMNS. C CALL SRTDAT(N,NNZ,INDROW,INDCOL,JPNTR,IWA) C C EXTEND THE DATA STRUCTURE TO ROWS. C CALL SETR(N,N,INDROW,JPNTR,INDCOL,IPNTR,IWA) C C DETERMINE THE DEGREE SEQUENCE FOR THE INTERSECTION C GRAPH OF THE COLUMNS OF L. C CALL DEGR(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(5*N+1),IWA(N+1)) C C COLOR THE INTERSECTION GRAPH OF THE COLUMNS OF L C WITH THE SMALLEST-LAST (SL) ORDERING. C CALL SLO(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(5*N+1),IWA(4*N+1), * MAXCLQ,IWA(1),IWA(N+1),IWA(2*N+1),IWA(3*N+1)) CALL SEQ(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(4*N+1),IWA(1), * MAXGRP,IWA(N+1)) DO 120 J = 1, N NGRP(J) = IWA(LISTP(J)) 120 CONTINUE C C EXIT IF THE SMALLEST-LAST ORDERING IS OPTIMAL. C IF (MAXGRP .EQ. MAXCLQ) GO TO 140 C C COLOR THE INTERSECTION GRAPH OF THE COLUMNS OF L C WITH THE INCIDENCE DEGREE (ID) ORDERING. C CALL IDO(N,N,INDROW,JPNTR,INDCOL,IPNTR,IWA(5*N+1),IWA(4*N+1), * MAXCLQ,IWA(1),IWA(N+1),IWA(2*N+1),IWA(3*N+1)) CALL SEQ(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(4*N+1),IWA(1), * NUMGRP,IWA(N+1)) C C RETAIN THE BETTER OF THE TWO ORDERINGS. C IF (NUMGRP .LT. MAXGRP) THEN MAXGRP = NUMGRP DO 130 J = 1, N NGRP(J) = IWA(LISTP(J)) 130 CONTINUE END IF 140 CONTINUE C C GENERATE THE SPARSITY PATTERN FOR THE LOWER C TRIANGULAR PART OF THE ORIGINAL MATRIX. C DO 150 J = 1, N IWA(LISTP(J)) = J 150 CONTINUE DO 170 J = 1, N DO 160 JP = JPNTR(J), JPNTR(J+1)-1 I = INDROW(JP) INDROW(JP) = MAX(IWA(I),IWA(J)) INDCOL(JP) = MIN(IWA(I),IWA(J)) 160 CONTINUE 170 CONTINUE C C SORT THE DATA STRUCTURE BY COLUMNS. C CALL SRTDAT(N,NNZ,INDROW,INDCOL,JPNTR,IWA) C C EXTEND THE DATA STRUCTURE TO ROWS. C CALL SETR(N,N,INDROW,JPNTR,INDCOL,IPNTR,IWA) END IF RETURN C C LAST CARD OF SUBROUTINE DSSM. C END SUBROUTINE IDOG(N,NGHBRP,NPNTRP,NGHBRS,NPNTRS,LISTP, * MAXCLQ,MAXID,IWA1,IWA2,IWA3) INTEGER N,MAXCLQ,MAXID INTEGER NGHBRP(*),NPNTRP(N+1),NGHBRS(*),NPNTRS(N+1),LISTP(N), * IWA1(0:N-1),IWA2(N),IWA3(N) C ********** C C SUBROUTINE IDOG C C GIVEN A LOOPLESS GRAPH G = (V,E), THIS SUBROUTINE DETERMINES C THE INCIDENCE DEGREE ORDERING OF THE VERTICES OF G. C C THE INCIDENCE DEGREE ORDERING IS DETERMINED RECURSIVELY BY C LETTING LIST(K), K = 1,...,N BE A VERTEX WITH MAXIMAL C INCIDENCE TO THE SUBGRAPH SPANNED BY THE ORDERED VERTICES. C AMONG ALL THE VERTICES OF MAXIMAL INCIDENCE, A VERTEX OF C MAXIMAL DEGREE IS CHOSEN. THIS SUBROUTINE DETERMINES THE C INVERSE OF THE INCIDENCE DEGREE ORDERING, THAT IS, AN ARRAY C LISTP SUCH THAT LISTP(LIST(K)) = K FOR K = 1,2,...,N. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE IDOG(N,NGHBRP,NPNTRP,NGHBRS,NPNTRS,LISTP, C MAXCLQ,MAXID,IWA1,IWA2,IWA3) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF VERTICES OF G. C C NGHBRP IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE C PREDECESSOR ADJACENCY LISTS FOR THE GRAPH G. C C NPNTRP IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE PREDECESSOR ADJACENCY C LISTS IN NGHBRP. THE VERTICES PRECEDING AND ADJACENT C TO VERTEX J ARE C C NGHBRP(K), K = NPNTRP(J),...,NPNTRP(J+1)-1. C C NOTE THAT NPNTRP(N+1)-1 IS THEN THE NUMBER OF VERTICES C PLUS EDGES OF THE GRAPH G. C C NGHBRS IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE C SUCCESSOR ADJACENCY LISTS FOR THE GRAPH G. C C NPNTRS IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE SUCCESSOR ADJACENCY C LISTS IN NGHBRS. THE VERTICES SUCCEEDING AND ADJACENT C TO VERTEX J ARE C C NGHBRS(K), K = NPNTRS(J),...,NPNTRS(J+1)-1. C C NOTE THAT NPNTRS(N+1)-1 IS THEN THE NUMBER OF VERTICES C PLUS EDGES OF THE GRAPH G. C C LISTP IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE INVERSE OF THE INCIDENCE DEGREE ORDERING OF THE C VERTICES. VERTEX J IS IN POSITION LISTP(J) OF THIS ORDERING. C C MAXCLQ IS AN INTEGER OUTPUT VARIABLE SET TO THE SIZE C OF THE LARGEST CLIQUE FOUND DURING THE ORDERING. C C MAXID IS AN INTEGER OUTPUT VARIABLE SET TO THE MAXIMUM C INCIDENCE DEGREE FOUND DURING THE ORDERING. C C IWA1,IWA2, AND IWA3 ARE INTEGER WORK ARRAYS OF LENGTH N. C C SUBPROGRAMS CALLED C C MINPACK-SUPPLIED ... NUMSRT C C FORTRAN-SUPPLIED ... MAX C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. DECEMBER 1984. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER I,J,K,MAXINC,MAXDEG,MAXLST,NCOMP,NUMDEG,NUMINC,NUMORD C C INITIALIZATION BLOCK. C DO 10 J = 1, N LISTP(J) = (NPNTRP(J+1) - NPNTRP(J) - 1) + * (NPNTRS(J+1) - NPNTRS(J) - 1) 10 CONTINUE MAXLST = (NPNTRP(N+1) + NPNTRS(N+1))/N C C SORT THE DEGREE SEQUENCE. C CALL NUMSRT(N,N-1,LISTP,1,IWA1,IWA2,IWA3) C C CREATE A DOUBLY-LINKED LIST TO ACCESS THE INCIDENCES OF THE C VERTICES. THE POINTERS FOR THE LINKED LIST ARE AS FOLLOWS. C C EACH UN-ORDERED VERTEX I IS IN A LIST (THE INCIDENCE LIST) C OF VERTICES WITH THE SAME INCIDENCE. C C IWA1(NUMINC) IS THE FIRST VERTEX IN THE NUMINC LIST C UNLESS IWA1(NUMINC) = 0. IN THIS CASE THERE ARE C NO VERTICES IN THE NUMINC LIST. C C IWA2(I) IS THE VERTEX BEFORE I IN THE INCIDENCE LIST C UNLESS IWA2(I) = 0. IN THIS CASE I IS THE FIRST C VERTEX IN THIS INCIDENCE LIST. C C IWA3(I) IS THE VERTEX AFTER I IN THE INCIDENCE LIST C UNLESS IWA3(I) = 0. IN THIS CASE I IS THE LAST C VERTEX IN THIS INCIDENCE LIST. C C IF I IS AN UN-ORDERED VERTEX, THEN -LISTP(I) IS THE C INCIDENCE OF I TO THE GRAPH INDUCED BY THE ORDERED C VERTICES. IF J IS AN ORDERED VERTEX, THEN LISTP(J) C IS THE INCIDENCE DEGREE ORDER OF VERTEX J. C MAXINC = 0 DO 20 J = 1, N I = IWA1(J-1) IWA1(J-1) = 0 IWA2(I) = 0 IWA3(I) = IWA1(0) IF (IWA1(0) .GT. 0) IWA2(IWA1(0)) = I IWA1(0) = I LISTP(J) = 0 20 CONTINUE MAXCLQ = 0 MAXID = 0 NUMORD = 1 C C BEGINNING OF ITERATION LOOP. C 30 CONTINUE C C CHOOSE A VERTEX J OF MAXIMAL DEGREE AMONG THE C VERTICES OF MAXIMAL INCIDENCE MAXINC. C 40 CONTINUE K = IWA1(MAXINC) IF (K .GT. 0) GO TO 50 MAXINC = MAXINC - 1 GO TO 40 50 CONTINUE MAXDEG = -1 DO 60 I = 1, MAXLST NUMDEG = (NPNTRP(K+1) - NPNTRP(K) - 1) + * (NPNTRS(K+1) - NPNTRS(K) - 1) IF (NUMDEG .GT. MAXDEG) THEN MAXDEG = NUMDEG J = K END IF K = IWA3(K) IF (K .LE. 0) GO TO 70 60 CONTINUE 70 CONTINUE LISTP(J) = NUMORD MAXID = MAX(MAXID,MAXINC) C C UPDATE THE SIZE OF THE LARGEST CLIQUE C FOUND DURING THE ORDERING. C IF (MAXINC .EQ. 0) NCOMP = 0 NCOMP = NCOMP + 1 IF (MAXINC+1 .EQ. NCOMP) MAXCLQ = MAX(MAXCLQ,NCOMP) C C TERMINATION TEST. C NUMORD = NUMORD + 1 IF (NUMORD .GT. N) GO TO 100 C C DELETE VERTEX J FROM THE MAXINC LIST. C IF (IWA2(J) .EQ. 0) THEN IWA1(MAXINC) = IWA3(J) ELSE IWA3(IWA2(J)) = IWA3(J) END IF IF (IWA3(J) .GT. 0) IWA2(IWA3(J)) = IWA2(J) C C DETERMINE ALL THE NEIGHBORS OF VERTEX J WHICH PRECEDE J C IN THE SUBGRAPH SPANNED BY THE UN-ORDERED VERTICES. C DO 80 K = NPNTRP(J), NPNTRP(J+1)-1 I = NGHBRP(K) C C UPDATE THE POINTERS TO THE CURRENT INCIDENCE LISTS. C NUMINC = -LISTP(I) IF (NUMINC .GE. 0) THEN LISTP(I) = LISTP(I) - 1 MAXINC = MAX(MAXINC,-LISTP(I)) C C DELETE VERTEX I FROM THE NUMINC LIST. C IF (IWA2(I) .EQ. 0) THEN IWA1(NUMINC) = IWA3(I) ELSE IWA3(IWA2(I)) = IWA3(I) END IF IF (IWA3(I) .GT. 0) IWA2(IWA3(I)) = IWA2(I) C C ADD VERTEX I TO THE NUMINC+1 LIST. C IWA2(I) = 0 IWA3(I) = IWA1(NUMINC+1) IF (IWA1(NUMINC+1) .GT. 0) IWA2(IWA1(NUMINC+1)) = I IWA1(NUMINC+1) = I END IF 80 CONTINUE C C DETERMINE ALL THE NEIGHBORS OF VERTEX J WHICH SUCCEED J C IN THE SUBGRAPH SPANNED BY THE UN-ORDERED VERTICES. C DO 90 K = NPNTRS(J), NPNTRS(J+1)-1 I = NGHBRS(K) C C UPDATE THE POINTERS TO THE CURRENT INCIDENCE LISTS. C NUMINC = -LISTP(I) IF (NUMINC .GE. 0) THEN LISTP(I) = LISTP(I) - 1 MAXINC = MAX(MAXINC,-LISTP(I)) C C DELETE VERTEX I FROM THE NUMINC LIST. C IF (IWA2(I) .EQ. 0) THEN IWA1(NUMINC) = IWA3(I) ELSE IWA3(IWA2(I)) = IWA3(I) END IF IF (IWA3(I) .GT. 0) IWA2(IWA3(I)) = IWA2(I) C C ADD VERTEX I TO THE NUMINC+1 LIST. C IWA2(I) = 0 IWA3(I) = IWA1(NUMINC+1) IF (IWA1(NUMINC+1) .GT. 0) IWA2(IWA1(NUMINC+1)) = I IWA1(NUMINC+1) = I END IF 90 CONTINUE C C END OF ITERATION LOOP. C GO TO 30 100 CONTINUE RETURN C C LAST CARD OF SUBROUTINE IDOG. C END SUBROUTINE SDPT(N,NGHBRP,NPNTRP,NGHBRS,NPNTRS,NGRP,MAXGRP, * IWA1,IWA2) INTEGER N,MAXGRP INTEGER NGHBRP(*),NPNTRP(N+1),NGHBRS(*),NPNTRS(N+1),NGRP(N), * IWA1(0:N-1),IWA2(N) C ********** C C SUBROUTINE SDPT C C GIVEN A LOOPLESS GRAPH G = (V,E), THIS SUBROUTINE DETERMINES C A SYMMETRIC COLORING OF G BY THE POWELL-TOINT DIRECT METHOD. C C THE POWELL-TOINT METHOD ASSIGNS THE K-TH COLOR BY EXAMINING C THE UN-COLORED VERTICES U(K) IN ORDER OF NON-INCREASING DEGREE C AND ASSIGNING COLOR K TO VERTEX V IF THERE ARE NO PATHS OF C LENGTH 1 OR 2 (IN THE GRAPH INDUCED BY U(K)) BETWEEN V AND C SOME K-COLORED VERTEX. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE SDPT(N,NGHBRP,NPNTRP,NGHBRS,NPNTRS,NGRP,MAXGRP, C IWA1,IWA2) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF VERTICES OF G. C C NGHBRP IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE C PREDECESSOR ADJACENCY LISTS FOR THE GRAPH G. C C NPNTRP IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE PREDECESSOR ADJACENCY C LISTS IN NGHBRP. THE VERTICES PRECEDING AND ADJACENT C TO VERTEX J ARE C C NGHBRP(K), K = NPNTRP(J),...,NPNTRP(J+1)-1. C C NOTE THAT NPNTRP(N+1)-1 IS THEN THE NUMBER OF VERTICES C PLUS EDGES OF THE GRAPH G. C C NGHBRS IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE C SUCCESSOR ADJACENCY LISTS FOR THE GRAPH G. C C NPNTRS IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE SUCCESSOR ADJACENCY C LISTS IN NGHBRS. THE VERTICES SUCCEEDING AND ADJACENT C TO VERTEX J ARE C C NGHBRS(K), K = NPNTRS(J),...,NPNTRS(J+1)-1. C C NOTE THAT NPNTRS(N+1)-1 IS THEN THE NUMBER OF VERTICES C PLUS EDGES OF THE GRAPH G. C C NGRP IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE SYMMETRIC COLORING OF G. VERTEX J IS COLORED WITH C COLOR NGRP(J). C C MAXGRP IS AN INTEGER OUTPUT VARIABLE WHICH SPECIFIES THE C NUMBER OF COLORS IN THE SYMMETRIC COLORING OF G. C C IWA1 AND IWA2 ARE INTEGER WORK ARRAYS OF LENGTH N. C C SUBPROGRAMS CALLED C C FORTRAN-SUPPLIED ... MAX C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. DECEMBER 1984. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER J,JP,K,KP,L,MAXDEG,NUMDEG,NUMV C C INITIALIZATION BLOCK. NUMV IS THE CURRENT NUMBER OF UN-COLORED C VERTICES, MAXDEG IS THE MAXIMUM INDUCED DEGREE OF THESE C VERTICES, AND MAXGRP IS THE CURRENT GROUP NUMBER (COLOR). C NUMV = N MAXDEG = 0 DO 10 J = 1, N NGRP(J) = (NPNTRP(J) - NPNTRP(J+1) + 1) + * (NPNTRS(J) - NPNTRS(J+1) + 1) MAXDEG = MAX(MAXDEG,-NGRP(J)) IWA2(J) = -J 10 CONTINUE MAXGRP = 0 C C BEGINNING OF ITERATION LOOP. C 20 CONTINUE C C SORT THE LIST OF UN-COLORED VERTICES SO THAT THEIR C INDUCED DEGREES ARE IN NON-DECREASING ORDER. C DO 30 NUMDEG = 0, MAXDEG IWA1(NUMDEG) = 0 30 CONTINUE DO 40 L = 1, NUMV NUMDEG = -NGRP(-IWA2(L)) IWA1(NUMDEG) = IWA1(NUMDEG) + 1 40 CONTINUE K = 1 DO 50 NUMDEG = MAXDEG, 0, -1 L = IWA1(NUMDEG) IWA1(NUMDEG) = K K = K + L 50 CONTINUE K = 1 60 CONTINUE J = IWA2(K) IF (J .GT. 0) THEN K = IWA1(-NGRP(J)) ELSE NUMDEG = -NGRP(-J) L = IWA1(NUMDEG) IWA2(K) = IWA2(L) IWA2(L) = -J IWA1(NUMDEG) = IWA1(NUMDEG) + 1 END IF IF (K .LE. NUMV) GO TO 60 MAXGRP = MAXGRP + 1 C C DETERMINE THE VERTICES IN GROUP MAXGRP. C DO 160 L = 1, NUMV J = IWA2(L) C C EXAMINE EACH VERTEX K PRECEDING VERTEX J AND ALL C THE NEIGHBORS OF VERTEX K TO DETERMINE IF VERTEX C J CAN BE CONSIDERED FOR GROUP MAXGRP. C DO 90 JP = NPNTRP(J), NPNTRP(J+1)-1 K = NGHBRP(JP) IF (NGRP(K) .EQ. MAXGRP) GO TO 150 IF (NGRP(K) .LE. 0) THEN DO 70 KP = NPNTRP(K), NPNTRP(K+1)-1 IF (NGRP(NGHBRP(KP)) .EQ. MAXGRP) GO TO 150 70 CONTINUE DO 80 KP = NPNTRS(K), NPNTRS(K+1)-1 IF (NGRP(NGHBRS(KP)) .EQ. MAXGRP) GO TO 150 80 CONTINUE END IF 90 CONTINUE C C EXAMINE EACH VERTEX K SUCCEEDING VERTEX J AND ALL C THE NEIGHBORS OF VERTEX K TO DETERMINE IF VERTEX C J CAN BE ADDED TO GROUP MAXGRP. C DO 120 JP = NPNTRS(J), NPNTRS(J+1)-1 K = NGHBRS(JP) IF (NGRP(K) .EQ. MAXGRP) GO TO 150 IF (NGRP(K) .LE. 0) THEN DO 100 KP = NPNTRP(K), NPNTRP(K+1)-1 IF (NGRP(NGHBRP(KP)) .EQ. MAXGRP) GO TO 150 100 CONTINUE DO 110 KP = NPNTRS(K), NPNTRS(K+1)-1 IF (NGRP(NGHBRS(KP)) .EQ. MAXGRP) GO TO 150 110 CONTINUE END IF 120 CONTINUE C C ADD VERTEX J TO GROUP MAXGRP AND REMOVE VERTEX J C FROM THE LIST OF UN-COLORED VERTICES. C NGRP(J) = MAXGRP IWA2(L) = 0 C C UPDATE THE DEGREES OF THE NEIGHBORS OF VERTEX J. C DO 130 JP = NPNTRP(J), NPNTRP(J+1)-1 K = NGHBRP(JP) IF (NGRP(K) .LT. 0) NGRP(K) = NGRP(K) + 1 130 CONTINUE DO 140 JP = NPNTRS(J), NPNTRS(J+1)-1 K = NGHBRS(JP) IF (NGRP(K) .LT. 0) NGRP(K) = NGRP(K) + 1 140 CONTINUE 150 CONTINUE 160 CONTINUE C C COMPRESS THE UPDATED LIST OF UN-COLORED VERTICES. C RESET NUMV AND RECOMPUTE MAXDEG. C K = 0 MAXDEG = 0 DO 170 L = 1, NUMV IF (IWA2(L) .NE. 0) THEN K = K + 1 IWA2(K) = -IWA2(L) MAXDEG = MAX(MAXDEG,-NGRP(IWA2(L))) END IF 170 CONTINUE NUMV = K C C END OF ITERATION LOOP. C IF (NUMV .GT. 0) GO TO 20 RETURN C C LAST CARD OF SUBROUTINE SDPT. C END SUBROUTINE SLOG(N,NGHBRP,NPNTRP,NGHBRS,NPNTRS,LISTP, * MAXCLQ,MAXVD,IWA1,IWA2,IWA3) INTEGER N,MAXCLQ,MAXVD INTEGER NGHBRP(*),NPNTRP(N+1),NGHBRS(*),NPNTRS(N+1),LISTP(N), * IWA1(0:N-1),IWA2(N),IWA3(N) C ********** C C SUBROUTINE SLOG C C GIVEN A LOOPLESS GRAPH G = (V,E), THIS SUBROUTINE DETERMINES C THE SMALLEST-LAST ORDERING OF THE VERTICES OF G. C C THE SMALLEST-LAST ORDERING IS DETERMINED RECURSIVELY BY C LETTING LIST(K), K = N,...,1 BE A VERTEX WITH LEAST DEGREE C IN THE SUBGRAPH SPANNED BY THE UN-ORDERED VERTICES. C THIS SUBROUTINE DETERMINES THE INVERSE OF THE SMALLEST-LAST C ORDERING, THAT IS, AN ARRAY LISTP SUCH THAT LISTP(LIST(K)) = K C FOR K = 1,2,...,N. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE SLOG(N,NGHBRP,NPNTRP,NGHBRS,NPNTRS,LISTP, C MAXCLQ,MAXVD,IWA1,IWA2,IWA3) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF VERTICES OF G. C C NGHBRP IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE C PREDECESSOR ADJACENCY LISTS FOR THE GRAPH G. C C NPNTRP IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE PREDECESSOR ADJACENCY C LISTS IN NGHBRP. THE VERTICES PRECEDING AND ADJACENT C TO VERTEX J ARE C C NGHBRP(K), K = NPNTRP(J),...,NPNTRP(J+1)-1. C C NOTE THAT NPNTRP(N+1)-1 IS THEN THE NUMBER OF VERTICES C PLUS EDGES OF THE GRAPH G. C C NGHBRS IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE C SUCCESSOR ADJACENCY LISTS FOR THE GRAPH G. C C NPNTRS IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE SUCCESSOR ADJACENCY C LISTS IN NGHBRS. THE VERTICES SUCCEEDING AND ADJACENT C TO VERTEX J ARE C C NGHBRS(K), K = NPNTRS(J),...,NPNTRS(J+1)-1. C C NOTE THAT NPNTRS(N+1)-1 IS THEN THE NUMBER OF VERTICES C PLUS EDGES OF THE GRAPH G. C C LISTP IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE INVERSE OF THE SMALLEST-LAST ORDERING OF THE VERTICES. C VERTEX J IS IN POSITION LISTP(J) OF THIS ORDERING. C C MAXCLQ IS AN INTEGER OUTPUT VARIABLE SET TO THE SIZE C OF THE LARGEST CLIQUE FOUND DURING THE ORDERING. C C MAXVD IS AN INTEGER OUTPUT VARIABLE SET TO THE MAXIMUM C VERTEX DEGREE FOUND DURING THE ORDERING. C C IWA1,IWA2, AND IWA3 ARE INTEGER WORK ARRAYS OF LENGTH N. C C SUBPROGRAMS CALLED C C FORTRAN-SUPPLIED ... MAX,MIN C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. DECEMBER 1984. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER I,J,K,MINDEG,NUMDEG,NUMORD C C INITIALIZATION BLOCK. C MINDEG = N DO 10 J = 1, N IWA1(J-1) = 0 LISTP(J) = (NPNTRP(J) - NPNTRP(J+1) + 1) + * (NPNTRS(J) - NPNTRS(J+1) + 1) MINDEG = MIN(MINDEG,-LISTP(J)) 10 CONTINUE C C CREATE A DOUBLY-LINKED LIST TO ACCESS THE DEGREES OF THE C VERTICES. THE POINTERS FOR THE LINKED LIST ARE AS FOLLOWS. C C EACH UN-ORDERED VERTEX I IS IN A LIST (THE DEGREE LIST) C OF VERTICES WITH THE SAME DEGREE. C C IWA1(NUMDEG) IS THE FIRST VERTEX IN THE NUMDEG LIST C UNLESS IWA1(NUMDEG) = 0. IN THIS CASE THERE ARE C NO VERTICES IN THE NUMDEG LIST. C C IWA2(I) IS THE VERTEX BEFORE I IN THE DEGREE LIST C UNLESS IWA2(I) = 0. IN THIS CASE I IS THE FIRST C VERTEX IN THIS DEGREE LIST. C C IWA3(I) IS THE VERTEX AFTER I IN THE DEGREE LIST C UNLESS IWA3(I) = 0. IN THIS CASE I IS THE LAST C VERTEX IN THIS DEGREE LIST. C C IF I IS AN UN-ORDERED VERTEX, THEN -LISTP(I) IS THE C DEGREE OF I IN THE GRAPH INDUCED BY THE UN-ORDERED C VERTICES. IF J IS AN ORDERED VERTEX, THEN LISTP(J) C IS THE SMALLEST-LAST ORDER OF VERTEX J. C DO 20 J = 1, N NUMDEG = -LISTP(J) IWA2(J) = 0 IWA3(J) = IWA1(NUMDEG) IF (IWA1(NUMDEG) .GT. 0) IWA2(IWA1(NUMDEG)) = J IWA1(NUMDEG) = J 20 CONTINUE MAXCLQ = 0 MAXVD = 0 NUMORD = N C C BEGINNING OF ITERATION LOOP. C 30 CONTINUE C C CHOOSE A VERTEX J OF MINIMAL DEGREE MINDEG. C 40 CONTINUE J = IWA1(MINDEG) IF (J .GT. 0) GO TO 50 MINDEG = MINDEG + 1 GO TO 40 50 CONTINUE LISTP(J) = NUMORD MAXVD = MAX(MAXVD,MINDEG) C C MARK THE SIZE OF THE LARGEST CLIQUE C FOUND DURING THE ORDERING. C IF (MINDEG+1 .EQ. NUMORD .AND. MAXCLQ .EQ. 0) * MAXCLQ = NUMORD C C TERMINATION TEST. C NUMORD = NUMORD - 1 IF (NUMORD .EQ. 0) GO TO 80 C C DELETE VERTEX J FROM THE MINDEG LIST. C IWA1(MINDEG) = IWA3(J) IF (IWA3(J) .GT. 0) IWA2(IWA3(J)) = 0 C C DETERMINE ALL THE NEIGHBORS OF VERTEX J WHICH PRECEDE J C IN THE SUBGRAPH SPANNED BY THE UN-ORDERED VERTICES. C DO 60 K = NPNTRP(J), NPNTRP(J+1)-1 I = NGHBRP(K) C C UPDATE THE POINTERS TO THE CURRENT DEGREE LISTS. C NUMDEG = -LISTP(I) IF (NUMDEG .GE. 0) THEN LISTP(I) = LISTP(I) + 1 MINDEG = MIN(MINDEG,-LISTP(I)) C C DELETE VERTEX I FROM THE NUMDEG LIST. C IF (IWA2(I) .EQ. 0) THEN IWA1(NUMDEG) = IWA3(I) ELSE IWA3(IWA2(I)) = IWA3(I) END IF IF (IWA3(I) .GT. 0) IWA2(IWA3(I)) = IWA2(I) C C ADD VERTEX I TO THE NUMDEG-1 LIST. C IWA2(I) = 0 IWA3(I) = IWA1(NUMDEG-1) IF (IWA1(NUMDEG-1) .GT. 0) IWA2(IWA1(NUMDEG-1)) = I IWA1(NUMDEG-1) = I END IF 60 CONTINUE C C DETERMINE ALL THE NEIGHBORS OF VERTEX J WHICH SUCCEED J C IN THE SUBGRAPH SPANNED BY THE UN-ORDERED VERTICES. C DO 70 K = NPNTRS(J), NPNTRS(J+1)-1 I = NGHBRS(K) C C UPDATE THE POINTERS TO THE CURRENT DEGREE LISTS. C NUMDEG = -LISTP(I) IF (NUMDEG .GE. 0) THEN LISTP(I) = LISTP(I) + 1 MINDEG = MIN(MINDEG,-LISTP(I)) C C DELETE VERTEX I FROM THE NUMDEG LIST. C IF (IWA2(I) .EQ. 0) THEN IWA1(NUMDEG) = IWA3(I) ELSE IWA3(IWA2(I)) = IWA3(I) END IF IF (IWA3(I) .GT. 0) IWA2(IWA3(I)) = IWA2(I) C C ADD VERTEX I TO THE NUMDEG-1 LIST. C IWA2(I) = 0 IWA3(I) = IWA1(NUMDEG-1) IF (IWA1(NUMDEG-1) .GT. 0) IWA2(IWA1(NUMDEG-1)) = I IWA1(NUMDEG-1) = I END IF 70 CONTINUE C C END OF ITERATION LOOP. C GO TO 30 80 CONTINUE RETURN C C LAST CARD OF SUBROUTINE SLOG. C END SUBROUTINE FDHS(N,INDROW,JPNTR,INDCOL,IPNTR,LISTP,NGRP, * MAXGRP,NUMGRP,ETA,FHESD,FHES,IWA) INTEGER N,MAXGRP,NUMGRP INTEGER INDROW(*),JPNTR(N+1),INDCOL(*),IPNTR(N+1), * LISTP(N),NGRP(N),IWA(N) REAL ETA(N),FHESD(N),FHES(*) C ********** C C SUBROUTINE FDHS C C THIS SUBROUTINE COMPUTES AN APPROXIMATION TO THE (SYMMETRIC) C HESSIAN MATRIX OF A FUNCTION BY A SUBSTITUTION METHOD. C THE LOWER TRIANGULAR PART OF THE APPROXIMATION IS STORED C WITH A COLUMN-ORIENTED DEFINITION OF THE SPARSITY PATTERN. C C THIS SUBROUTINE REQUIRES A SYMMETRIC PERMUTATION OF THE C HESSIAN MATRIX AND A PARTITION OF THE COLUMNS OF THE HESSIAN C MATRIX CONSISTENT WITH THE DETERMINATION OF THE HESSIAN C MATRIX BY A LOWER TRIANGULAR SUBSTITUTION METHOD. C THIS INFORMATION CAN BE PROVIDED BY SUBROUTINE DSSM. C C THE SYMMETRIC PERMUTATION OF THE HESSIAN MATRIX IS DEFINED C BY THE ARRAY LISTP. THIS ARRAY IS ONLY USED INTERNALLY. C C THE PARTITION OF THE HESSIAN MATRIX IS DEFINED BY THE ARRAY C NGRP BY SETTING NGRP(J) TO THE GROUP NUMBER OF COLUMN J. C THE USER MUST PROVIDE AN APPROXIMATION TO THE COLUMNS OF C THE HESSIAN MATRIX IN EACH GROUP BY SPECIFYING A DIFFERENCE C PARAMETER VECTOR ETA AND AN APPROXIMATION TO H*D WHERE H IS C THE HESSIAN MATRIX AND THE VECTOR D IS DEFINED BY THE C FOLLOWING SECTION OF CODE. C C DO 10 J = 1, N C D(J) = 0.0 C IF (NGRP(J) .EQ. NUMGRP) D(J) = ETA(J) C 10 CONTINUE C C IN THE ABOVE CODE NUMGRP IS A GROUP NUMBER AND ETA(J) IS THE C DIFFERENCE PARAMETER USED TO APPROXIMATE COLUMN J OF THE C HESSIAN MATRIX. SUITABLE VALUES FOR ETA(J) MUST BE PROVIDED. C C AS MENTIONED ABOVE, AN APPROXIMATION TO H*D MUST BE PROVIDED. C FOR EXAMPLE, IF GRAD(X) IS THE GRADIENT OF THE FUNCTION AT X, C THEN C C GRAD(X+D) - GRAD(X) C C CORRESPONDS TO THE FORWARD DIFFERENCE APPROXIMATION. C C THE LOWER TRIANGULAR SUBSTITUTION METHOD REQUIRES THAT THE C APPROXIMATIONS TO H*D FOR ALL THE GROUPS BE STORED IN SPECIAL C LOCATIONS OF THE ARRAY FHES. THIS IS DONE BY CALLING FDHS C SUCCESSIVELY WITH NUMGRP = 1,2,...,MAXGRP. ON THE CALL WITH C NUMGRP = MAXGRP, THE SUBROUTINE THEN PROCEEDS TO OVERWRITE C FHES WITH THE APPROXIMATION TO THE LOWER TRIANGULAR PART OF C THE HESSIAN MATRIX. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE FDHS(N,INDROW,JPNTR,INDCOL,IPNTR,LISTP,NGRP, C MAXGRP,NUMGRP,ETA,FHESD,FHES,IWA) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE ORDER C OF THE HESSIAN MATRIX. C C INDROW IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW C INDICES FOR THE NON-ZEROES IN THE LOWER TRIANGULAR PART C OF THE HESSIAN MATRIX. C C JPNTR IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW. C THE ROW INDICES FOR COLUMN J ARE C C INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1. C C NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZEROES C IN THE LOWER TRIANGULAR PART OF THE HESSIAN MATRIX. C C INDCOL IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE COLUMN C INDICES FOR THE NON-ZEROES IN THE LOWER TRIANGULAR PART C OF THE HESSIAN MATRIX. C C IPNTR IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL. C THE COLUMN INDICES FOR ROW I ARE C C INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1. C C NOTE THAT IPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZEROES C IN THE LOWER TRIANGULAR PART OF THE HESSIAN MATRIX. C C LISTP IS AN INTEGER INPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE SYMMETRIC PERMUTATION OF THE HESSIAN MATRIX. ELEMENT C (I,J) OF THE HESSIAN MATRIX IS THE (LISTP(I),LISTP(J)) C ELEMENT OF THE PERMUTED HESSIAN. C C NGRP IS AN INTEGER INPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE PARTITION OF THE COLUMNS OF THE HESSIAN MATRIX. C COLUMN J BELONGS TO GROUP NGRP(J). C C MAXGRP IS A POSITIVE INTEGER INPUT VARIABLE WHICH SPECIFIES C THE NUMBER OF GROUPS IN THE PARTITION OF THE COLUMNS OF C THE HESSIAN MATRIX. C C NUMGRP IS A POSITIVE INTEGER INPUT VARIABLE SET TO A GROUP C NUMBER IN THE PARTITION. C C ETA IS AN INPUT ARRAY OF LENGTH N WHICH CONTAINS THE C DIFFERENCE PARAMETER VECTOR. C C FHESD IS AN INPUT ARRAY OF LENGTH N WHICH CONTAINS AN C APPROXIMATION TO H*D, WHERE H IS THE HESSIAN MATRIX C AND D IS THE DIFFERENCE VECTOR FOR GROUP NUMGRP. C C FHES IS AN OUTPUT ARRAY OF LENGTH NNZ, WHERE NNZ IS THE C NUMBER OF NON-ZERO ELEMENTS IN THE LOWER TRIANGULAR PART C OF THE HESSIAN MATRIX. ON OUTPUT WITH NUMGRP LESS THAN C MAXGRP, THE FHESD ARRAY FOR GROUP NUMGRP HAS BEEN STORED C IN FHES. WHEN NUMGRP = MAXGRP THE SUBROUTINE OVERWRITES C FHES WITH AN APPROXIMATION TO THE LOWER TRIANGULAR PART C OF THE HESSIAN MATRIX. THE APPROXIMATION IS STORED IN C FHES WITH A COLUMN-ORIENTED DEFINITION OF THE SPARSITY C PATTERN. THUS THE ELEMENTS IN COLUMN J OF THE LOWER C TRIANGULAR PART OF THE HESSIAN MATRIX ARE C C FHES(K), K = JPNTR(J),...,JPNTR(J+1)-1, C C AND THE ROW INDICES FOR THESE ELEMENTS ARE C C INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1. C C IWA IS AN INTEGER WORK ARRAY OF LENGTH N. C C SUBPROGRAMS CALLED C C FORTRAN-SUPPLIED ... ABS C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. DECEMBER 1984. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER I,IP,IROW,J,JP,K,L,NUMG,NUML REAL SUM C C STORE THE I-TH ELEMENT OF GRADIENT DIFFERENCE FHESD C CORRESPONDING TO GROUP NUMGRP IF THERE IS A POSITION C (I,J) SUCH THAT NGRP(J) = NUMGRP AND (I,J) IS MAPPED C ONTO THE LOWER TRIANGULAR PART OF THE PERMUTED MATRIX. C DO 50 J = 1, N IF (NGRP(J) .EQ. NUMGRP) THEN NUML = LISTP(J) DO 30 IP = IPNTR(J), IPNTR(J+1)-1 I = INDCOL(IP) IF (LISTP(I) .GT. NUML) THEN DO 10 JP = JPNTR(I), JPNTR(I+1)-1 IF (INDROW(JP) .EQ. J) THEN FHES(JP) = FHESD(I) GO TO 20 END IF 10 CONTINUE 20 CONTINUE END IF 30 CONTINUE DO 40 JP = JPNTR(J), JPNTR(J+1)-1 I = INDROW(JP) IF (LISTP(I) .GE. NUML) FHES(JP) = FHESD(I) 40 CONTINUE END IF 50 CONTINUE C C EXIT IF THIS IS NOT THE LAST GROUP. C IF (NUMGRP .LT. MAXGRP) RETURN C C MARK ALL COLUMN INDICES J SUCH THAT (I,J) IS MAPPED ONTO C THE LOWER TRIANGULAR PART OF THE PERMUTED MATRIX. C DO 80 I = 1, N NUML = LISTP(I) DO 60 IP = IPNTR(I), IPNTR(I+1)-1 J = INDCOL(IP) IF (NUML .GE. LISTP(J)) INDCOL(IP) = -INDCOL(IP) 60 CONTINUE DO 70 JP = JPNTR(I), JPNTR(I+1)-1 J = INDROW(JP) IF (NUML .GT. LISTP(J)) INDROW(JP) = -INDROW(JP) 70 CONTINUE 80 CONTINUE C C INVERT THE ARRAY LISTP. C DO 90 J = 1, N IWA(LISTP(J)) = J 90 CONTINUE DO 100 J = 1, N LISTP(J) = IWA(J) 100 CONTINUE C C DETERMINE THE LOWER TRIANGULAR PART OF THE ORIGINAL MATRIX. C DO 220 IROW = N, 1, -1 I = LISTP(IROW) C C FIND THE POSITIONS OF THE ELEMENTS IN THE I-TH ROW OF THE C LOWER TRIANGULAR PART OF THE ORIGINAL MATRIX THAT HAVE C ALREADY BEEN DETERMINED. C DO 130 IP = IPNTR(I), IPNTR(I+1)-1 J = INDCOL(IP) IF (J .GT. 0) THEN DO 110 JP = JPNTR(J), JPNTR(J+1)-1 IF (INDROW(JP) .EQ. I) THEN IWA(J) = JP GO TO 120 END IF 110 CONTINUE 120 CONTINUE END IF 130 CONTINUE C C DETERMINE THE ELEMENTS IN THE I-TH ROW OF THE LOWER C TRIANGULAR PART OF THE ORIGINAL MATRIX WHICH GET MAPPED C ONTO THE LOWER TRIANGULAR PART OF THE PERMUTED MATRIX. C DO 180 K = IPNTR(I), IPNTR(I+1)-1 J = -INDCOL(K) IF (J .GT. 0) THEN INDCOL(K) = J C C DETERMINE THE (I,J) ELEMENT. C NUMG = NGRP(J) SUM = 0.0 DO 140 IP = IPNTR(I), IPNTR(I+1)-1 L = ABS(INDCOL(IP)) IF (NGRP(L) .EQ. NUMG .AND. L .NE. J) * SUM = SUM + FHES(IWA(L))*ETA(L) 140 CONTINUE DO 150 JP = JPNTR(I), JPNTR(I+1)-1 L = ABS(INDROW(JP)) IF (NGRP(L) .EQ. NUMG .AND. L .NE. J) * SUM = SUM + FHES(JP)*ETA(L) 150 CONTINUE C C STORE THE (I,J) ELEMENT. C DO 160 JP = JPNTR(J), JPNTR(J+1)-1 IF (INDROW(JP) .EQ. I) THEN FHES(JP) = (FHES(JP) - SUM)/ETA(J) GO TO 170 END IF 160 CONTINUE 170 CONTINUE END IF 180 CONTINUE C C DETERMINE THE ELEMENTS IN THE I-TH ROW OF THE STRICT UPPER C TRIANGULAR PART OF THE ORIGINAL MATRIX WHICH GET MAPPED C ONTO THE LOWER TRIANGULAR PART OF THE PERMUTED MATRIX. C DO 210 K = JPNTR(I), JPNTR(I+1)-1 J = -INDROW(K) IF (J .GT. 0) THEN INDROW(K) = J C C DETERMINE THE (I,J) ELEMENT. C NUMG = NGRP(J) SUM = 0.0 DO 190 IP = IPNTR(I), IPNTR(I+1)-1 L = ABS(INDCOL(IP)) IF (NGRP(L) .EQ. NUMG) * SUM = SUM + FHES(IWA(L))*ETA(L) 190 CONTINUE DO 200 JP = JPNTR(I), JPNTR(I+1)-1 L = ABS(INDROW(JP)) IF (NGRP(L) .EQ. NUMG .AND. L .NE. J) * SUM = SUM + FHES(JP)*ETA(L) 200 CONTINUE C C STORE THE (I,J) ELEMENT. C FHES(K) = (FHES(K) - SUM)/ETA(J) END IF 210 CONTINUE 220 CONTINUE C C RE-INVERT THE ARRAY LISTP. C DO 230 J = 1, N IWA(LISTP(J)) = J 230 CONTINUE DO 240 J = 1, N LISTP(J) = IWA(J) 240 CONTINUE RETURN C C LAST CARD OF SUBROUTINE FDHS. C END SUBROUTINE DEGR(N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,IWA) INTEGER N INTEGER INDROW(*),JPNTR(N+1),INDCOL(*),IPNTR(*),NDEG(N),IWA(N) C ********** C C SUBROUTINE DEGR C C GIVEN THE SPARSITY PATTERN OF AN M BY N MATRIX A, C THIS SUBROUTINE DETERMINES THE DEGREE SEQUENCE FOR C THE INTERSECTION GRAPH OF THE COLUMNS OF A. C C IN GRAPH-THEORY TERMINOLOGY, THE INTERSECTION GRAPH OF C THE COLUMNS OF A IS THE LOOPLESS GRAPH G WITH VERTICES C A(J), J = 1,2,...,N WHERE A(J) IS THE J-TH COLUMN OF A C AND WITH EDGE (A(I),A(J)) IF AND ONLY IF COLUMNS I AND J C HAVE A NON-ZERO IN THE SAME ROW POSITION. C C NOTE THAT THE VALUE OF M IS NOT NEEDED BY DEGR AND IS C THEREFORE NOT PRESENT IN THE SUBROUTINE STATEMENT. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE DEGR(N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,IWA) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF COLUMNS OF A. C C INDROW IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW C INDICES FOR THE NON-ZEROES IN THE MATRIX A. C C JPNTR IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW. C THE ROW INDICES FOR COLUMN J ARE C C INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1. C C NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO C ELEMENTS OF THE MATRIX A. C C INDCOL IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE C COLUMN INDICES FOR THE NON-ZEROES IN THE MATRIX A. C C IPNTR IS AN INTEGER INPUT ARRAY OF LENGTH M + 1 WHICH C SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL. C THE COLUMN INDICES FOR ROW I ARE C C INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1. C C NOTE THAT IPNTR(M+1)-1 IS THEN THE NUMBER OF NON-ZERO C ELEMENTS OF THE MATRIX A. C C NDEG IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH C SPECIFIES THE DEGREE SEQUENCE. THE DEGREE OF THE C J-TH COLUMN OF A IS NDEG(J). C C IWA IS AN INTEGER WORK ARRAY OF LENGTH N. C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER IC,IP,IR,JCOL,JP C C INITIALIZATION BLOCK. C DO 10 JP = 1, N NDEG(JP) = 0 IWA(JP) = 0 10 CONTINUE C C COMPUTE THE DEGREE SEQUENCE BY DETERMINING THE CONTRIBUTIONS C TO THE DEGREES FROM THE CURRENT(JCOL) COLUMN AND FURTHER C COLUMNS WHICH HAVE NOT YET BEEN CONSIDERED. C DO 40 JCOL = 2, N IWA(JCOL) = N C C DETERMINE ALL POSITIONS (IR,JCOL) WHICH CORRESPOND C TO NON-ZEROES IN THE MATRIX. C DO 30 JP = JPNTR(JCOL), JPNTR(JCOL+1)-1 IR = INDROW(JP) C C FOR EACH ROW IR, DETERMINE ALL POSITIONS (IR,IC) C WHICH CORRESPOND TO NON-ZEROES IN THE MATRIX. C DO 20 IP = IPNTR(IR), IPNTR(IR+1)-1 IC = INDCOL(IP) C C ARRAY IWA MARKS COLUMNS WHICH HAVE CONTRIBUTED TO C THE DEGREE COUNT OF COLUMN JCOL. UPDATE THE DEGREE C COUNTS OF THESE COLUMNS AS WELL AS COLUMN JCOL. C IF (IWA(IC) .LT. JCOL) THEN IWA(IC) = JCOL NDEG(IC) = NDEG(IC) + 1 NDEG(JCOL) = NDEG(JCOL) + 1 END IF 20 CONTINUE 30 CONTINUE 40 CONTINUE RETURN C C LAST CARD OF SUBROUTINE DEGR. C END SUBROUTINE IDO(M,N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,LIST, * MAXCLQ,IWA1,IWA2,IWA3,IWA4) INTEGER M,N,MAXCLQ INTEGER INDROW(*),JPNTR(N+1),INDCOL(*),IPNTR(M+1),NDEG(N), * LIST(N),IWA1(0:N-1),IWA2(N),IWA3(N),IWA4(N) C ********** C C SUBROUTINE IDO C C GIVEN THE SPARSITY PATTERN OF AN M BY N MATRIX A, THIS C SUBROUTINE DETERMINES AN INCIDENCE-DEGREE ORDERING OF THE C COLUMNS OF A. C C THE INCIDENCE-DEGREE ORDERING IS DEFINED FOR THE LOOPLESS C GRAPH G WITH VERTICES A(J), J = 1,2,...,N WHERE A(J) IS THE C J-TH COLUMN OF A AND WITH EDGE (A(I),A(J)) IF AND ONLY IF C COLUMNS I AND J HAVE A NON-ZERO IN THE SAME ROW POSITION. C C THE INCIDENCE-DEGREE ORDERING IS DETERMINED RECURSIVELY BY C LETTING LIST(K), K = 1,...,N BE A COLUMN WITH MAXIMAL C INCIDENCE TO THE SUBGRAPH SPANNED BY THE ORDERED COLUMNS. C AMONG ALL THE COLUMNS OF MAXIMAL INCIDENCE, IDO CHOOSES A C COLUMN OF MAXIMAL DEGREE. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE IDO(M,N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,LIST, C MAXCLQ,IWA1,IWA2,IWA3,IWA4) C C WHERE C C M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF ROWS OF A. C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF COLUMNS OF A. C C INDROW IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW C INDICES FOR THE NON-ZEROES IN THE MATRIX A. C C JPNTR IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW. C THE ROW INDICES FOR COLUMN J ARE C C INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1. C C NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO C ELEMENTS OF THE MATRIX A. C C INDCOL IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE C COLUMN INDICES FOR THE NON-ZEROES IN THE MATRIX A. C C IPNTR IS AN INTEGER INPUT ARRAY OF LENGTH M + 1 WHICH C SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL. C THE COLUMN INDICES FOR ROW I ARE C C INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1. C C NOTE THAT IPNTR(M+1)-1 IS THEN THE NUMBER OF NON-ZERO C ELEMENTS OF THE MATRIX A. C C NDEG IS AN INTEGER INPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE DEGREE SEQUENCE. THE DEGREE OF THE J-TH COLUMN C OF A IS NDEG(J). C C LIST IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE INCIDENCE-DEGREE ORDERING OF THE COLUMNS OF A. THE J-TH C COLUMN IN THIS ORDER IS LIST(J). C C MAXCLQ IS AN INTEGER OUTPUT VARIABLE SET TO THE SIZE C OF THE LARGEST CLIQUE FOUND DURING THE ORDERING. C C IWA1,IWA2,IWA3, AND IWA4 ARE INTEGER WORK ARRAYS OF LENGTH N. C C SUBPROGRAMS CALLED C C MINPACK-SUPPLIED ... NUMSRT C C FORTRAN-SUPPLIED ... MAX C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. AUGUST 1984. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER IC,IP,IR,JCOL,JP, * MAXINC,MAXLST,NCOMP,NUMINC,NUMLST,NUMORD,NUMWGT C C SORT THE DEGREE SEQUENCE. C CALL NUMSRT(N,N-1,NDEG,-1,IWA4,IWA2,IWA3) C C INITIALIZATION BLOCK. C C CREATE A DOUBLY-LINKED LIST TO ACCESS THE INCIDENCES OF THE C COLUMNS. THE POINTERS FOR THE LINKED LIST ARE AS FOLLOWS. C C EACH UN-ORDERED COLUMN IC IS IN A LIST (THE INCIDENCE LIST) C OF COLUMNS WITH THE SAME INCIDENCE. C C IWA1(NUMINC) IS THE FIRST COLUMN IN THE NUMINC LIST C UNLESS IWA1(NUMINC) = 0. IN THIS CASE THERE ARE C NO COLUMNS IN THE NUMINC LIST. C C IWA2(IC) IS THE COLUMN BEFORE IC IN THE INCIDENCE LIST C UNLESS IWA2(IC) = 0. IN THIS CASE IC IS THE FIRST C COLUMN IN THIS INCIDENCE LIST. C C IWA3(IC) IS THE COLUMN AFTER IC IN THE INCIDENCE LIST C UNLESS IWA3(IC) = 0. IN THIS CASE IC IS THE LAST C COLUMN IN THIS INCIDENCE LIST. C C IF IC IS AN UN-ORDERED COLUMN, THEN LIST(IC) IS THE C INCIDENCE OF IC TO THE GRAPH INDUCED BY THE ORDERED C COLUMNS. IF JCOL IS AN ORDERED COLUMN, THEN LIST(JCOL) C IS THE INCIDENCE-DEGREE ORDER OF COLUMN JCOL. C MAXINC = 0 DO 10 JP = N, 1, -1 IC = IWA4(JP) IWA1(N-JP) = 0 IWA2(IC) = 0 IWA3(IC) = IWA1(0) IF (IWA1(0) .GT. 0) IWA2(IWA1(0)) = IC IWA1(0) = IC IWA4(JP) = 0 LIST(JP) = 0 10 CONTINUE C C DETERMINE THE MAXIMAL SEARCH LENGTH FOR THE LIST C OF COLUMNS OF MAXIMAL INCIDENCE. C MAXLST = 0 DO 20 IR = 1, M MAXLST = MAXLST + (IPNTR(IR+1) - IPNTR(IR))**2 20 CONTINUE MAXLST = MAXLST/N MAXCLQ = 0 NUMORD = 1 C C BEGINNING OF ITERATION LOOP. C 30 CONTINUE C C CHOOSE A COLUMN JCOL OF MAXIMAL DEGREE AMONG THE C COLUMNS OF MAXIMAL INCIDENCE MAXINC. C 40 CONTINUE JP = IWA1(MAXINC) IF (JP .GT. 0) GO TO 50 MAXINC = MAXINC - 1 GO TO 40 50 CONTINUE NUMWGT = -1 DO 60 NUMLST = 1, MAXLST IF (NDEG(JP) .GT. NUMWGT) THEN NUMWGT = NDEG(JP) JCOL = JP END IF JP = IWA3(JP) IF (JP .LE. 0) GO TO 70 60 CONTINUE 70 CONTINUE LIST(JCOL) = NUMORD C C UPDATE THE SIZE OF THE LARGEST CLIQUE C FOUND DURING THE ORDERING. C IF (MAXINC .EQ. 0) NCOMP = 0 NCOMP = NCOMP + 1 IF (MAXINC+1 .EQ. NCOMP) MAXCLQ = MAX(MAXCLQ,NCOMP) C C TERMINATION TEST. C NUMORD = NUMORD + 1 IF (NUMORD .GT. N) GO TO 100 C C DELETE COLUMN JCOL FROM THE MAXINC LIST. C IF (IWA2(JCOL) .EQ. 0) THEN IWA1(MAXINC) = IWA3(JCOL) ELSE IWA3(IWA2(JCOL)) = IWA3(JCOL) END IF IF (IWA3(JCOL) .GT. 0) IWA2(IWA3(JCOL)) = IWA2(JCOL) C C FIND ALL COLUMNS ADJACENT TO COLUMN JCOL. C IWA4(JCOL) = N C C DETERMINE ALL POSITIONS (IR,JCOL) WHICH CORRESPOND C TO NON-ZEROES IN THE MATRIX. C DO 90 JP = JPNTR(JCOL), JPNTR(JCOL+1)-1 IR = INDROW(JP) C C FOR EACH ROW IR, DETERMINE ALL POSITIONS (IR,IC) C WHICH CORRESPOND TO NON-ZEROES IN THE MATRIX. C DO 80 IP = IPNTR(IR), IPNTR(IR+1)-1 IC = INDCOL(IP) C C ARRAY IWA4 MARKS COLUMNS WHICH ARE ADJACENT TO C COLUMN JCOL. C IF (IWA4(IC) .LT. NUMORD) THEN IWA4(IC) = NUMORD C C UPDATE THE POINTERS TO THE CURRENT INCIDENCE LISTS. C NUMINC = LIST(IC) LIST(IC) = LIST(IC) + 1 MAXINC = MAX(MAXINC,LIST(IC)) C C DELETE COLUMN IC FROM THE NUMINC LIST. C IF (IWA2(IC) .EQ. 0) THEN IWA1(NUMINC) = IWA3(IC) ELSE IWA3(IWA2(IC)) = IWA3(IC) END IF IF (IWA3(IC) .GT. 0) IWA2(IWA3(IC)) = IWA2(IC) C C ADD COLUMN IC TO THE NUMINC+1 LIST. C IWA2(IC) = 0 IWA3(IC) = IWA1(NUMINC+1) IF (IWA1(NUMINC+1) .GT. 0) IWA2(IWA1(NUMINC+1)) = IC IWA1(NUMINC+1) = IC END IF 80 CONTINUE 90 CONTINUE C C END OF ITERATION LOOP. C GO TO 30 100 CONTINUE C C INVERT THE ARRAY LIST. C DO 110 JCOL = 1, N IWA2(LIST(JCOL)) = JCOL 110 CONTINUE DO 120 JP = 1, N LIST(JP) = IWA2(JP) 120 CONTINUE RETURN C C LAST CARD OF SUBROUTINE IDO. C END SUBROUTINE NUMSRT(N,NMAX,NUM,MODE,INDEX,LAST,NEXT) INTEGER N,NMAX,MODE INTEGER NUM(N),INDEX(N),LAST(0:NMAX),NEXT(N) C **********. C C SUBROUTINE NUMSRT C C GIVEN A SEQUENCE OF INTEGERS, THIS SUBROUTINE GROUPS C TOGETHER THOSE INDICES WITH THE SAME SEQUENCE VALUE C AND, OPTIONALLY, SORTS THE SEQUENCE INTO EITHER C ASCENDING OR DESCENDING ORDER. C C THE SEQUENCE OF INTEGERS IS DEFINED BY THE ARRAY NUM, C AND IT IS ASSUMED THAT THE INTEGERS ARE EACH FROM THE SET C 0,1,...,NMAX. ON OUTPUT THE INDICES K SUCH THAT NUM(K) = L C FOR ANY L = 0,1,...,NMAX CAN BE OBTAINED FROM THE ARRAYS C LAST AND NEXT AS FOLLOWS. C C K = LAST(L) C WHILE (K .NE. 0) K = NEXT(K) C C OPTIONALLY, THE SUBROUTINE PRODUCES AN ARRAY INDEX SO THAT C THE SEQUENCE NUM(INDEX(I)), I = 1,2,...,N IS SORTED. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE NUMSRT(N,NMAX,NUM,MODE,INDEX,LAST,NEXT) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE. C C NMAX IS A POSITIVE INTEGER INPUT VARIABLE. C C NUM IS AN INPUT ARRAY OF LENGTH N WHICH CONTAINS THE C SEQUENCE OF INTEGERS TO BE GROUPED AND SORTED. IT C IS ASSUMED THAT THE INTEGERS ARE EACH FROM THE SET C 0,1,...,NMAX. C C MODE IS AN INTEGER INPUT VARIABLE. THE SEQUENCE NUM IS C SORTED IN ASCENDING ORDER IF MODE IS POSITIVE AND IN C DESCENDING ORDER IF MODE IS NEGATIVE. IF MODE IS 0, C NO SORTING IS DONE. C C INDEX IS AN INTEGER OUTPUT ARRAY OF LENGTH N SET SO C THAT THE SEQUENCE C C NUM(INDEX(I)), I = 1,2,...,N C C IS SORTED ACCORDING TO THE SETTING OF MODE. IF MODE C IS 0, INDEX IS NOT REFERENCED. C C LAST IS AN INTEGER OUTPUT ARRAY OF LENGTH NMAX + 1. THE C INDEX OF NUM FOR THE LAST OCCURRENCE OF L IS LAST(L) C FOR ANY L = 0,1,...,NMAX UNLESS LAST(L) = 0. IN C THIS CASE L DOES NOT APPEAR IN NUM. C C NEXT IS AN INTEGER OUTPUT ARRAY OF LENGTH N. IF C NUM(K) = L, THEN THE INDEX OF NUM FOR THE PREVIOUS C OCCURRENCE OF L IS NEXT(K) FOR ANY L = 0,1,...,NMAX C UNLESS NEXT(K) = 0. IN THIS CASE THERE IS NO PREVIOUS C OCCURRENCE OF L IN NUM. C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER I,J,JINC,JL,JU,K,L C C DETERMINE THE ARRAYS NEXT AND LAST. C DO 10 I = 0, NMAX LAST(I) = 0 10 CONTINUE DO 20 K = 1, N L = NUM(K) NEXT(K) = LAST(L) LAST(L) = K 20 CONTINUE IF (MODE .EQ. 0) RETURN C C STORE THE POINTERS TO THE SORTED ARRAY IN INDEX. C I = 1 IF (MODE .GT. 0) THEN JL = 0 JU = NMAX JINC = 1 ELSE JL = NMAX JU = 0 JINC = -1 END IF DO 50 J = JL, JU, JINC K = LAST(J) 30 CONTINUE IF (K .EQ. 0) GO TO 40 INDEX(I) = K I = I + 1 K = NEXT(K) GO TO 30 40 CONTINUE 50 CONTINUE RETURN C C LAST CARD OF SUBROUTINE NUMSRT. C END SUBROUTINE SEQ(N,INDROW,JPNTR,INDCOL,IPNTR,LIST,NGRP,MAXGRP, * IWA) INTEGER N,MAXGRP INTEGER INDROW(*),JPNTR(N+1),INDCOL(*),IPNTR(*),LIST(N), * NGRP(N),IWA(N) C ********** C C SUBROUTINE SEQ C C GIVEN THE SPARSITY PATTERN OF AN M BY N MATRIX A, THIS C SUBROUTINE DETERMINES A CONSISTENT PARTITION OF THE C COLUMNS OF A BY A SEQUENTIAL ALGORITHM. C C A CONSISTENT PARTITION IS DEFINED IN TERMS OF THE LOOPLESS C GRAPH G WITH VERTICES A(J), J = 1,2,...,N WHERE A(J) IS THE C J-TH COLUMN OF A AND WITH EDGE (A(I),A(J)) IF AND ONLY IF C COLUMNS I AND J HAVE A NON-ZERO IN THE SAME ROW POSITION. C C A PARTITION OF THE COLUMNS OF A INTO GROUPS IS CONSISTENT C IF THE COLUMNS IN ANY GROUP ARE NOT ADJACENT IN THE GRAPH G. C IN GRAPH-THEORY TERMINOLOGY, A CONSISTENT PARTITION OF THE C COLUMNS OF A CORRESPONDS TO A COLORING OF THE GRAPH G. C C THE SUBROUTINE EXAMINES THE COLUMNS IN THE ORDER SPECIFIED C BY THE ARRAY LIST, AND ASSIGNS THE CURRENT COLUMN TO THE C GROUP WITH THE SMALLEST POSSIBLE NUMBER. C C NOTE THAT THE VALUE OF M IS NOT NEEDED BY SEQ AND IS C THEREFORE NOT PRESENT IN THE SUBROUTINE STATEMENT. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE SEQ(N,INDROW,JPNTR,INDCOL,IPNTR,LIST,NGRP,MAXGRP, C IWA) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF COLUMNS OF A. C C INDROW IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW C INDICES FOR THE NON-ZEROES IN THE MATRIX A. C C JPNTR IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW. C THE ROW INDICES FOR COLUMN J ARE C C INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1. C C NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO C ELEMENTS OF THE MATRIX A. C C INDCOL IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE C COLUMN INDICES FOR THE NON-ZEROES IN THE MATRIX A. C C IPNTR IS AN INTEGER INPUT ARRAY OF LENGTH M + 1 WHICH C SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL. C THE COLUMN INDICES FOR ROW I ARE C C INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1. C C NOTE THAT IPNTR(M+1)-1 IS THEN THE NUMBER OF NON-ZERO C ELEMENTS OF THE MATRIX A. C C LIST IS AN INTEGER INPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE ORDER TO BE USED BY THE SEQUENTIAL ALGORITHM. C THE J-TH COLUMN IN THIS ORDER IS LIST(J). C C NGRP IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE PARTITION OF THE COLUMNS OF A. COLUMN JCOL BELONGS C TO GROUP NGRP(JCOL). C C MAXGRP IS AN INTEGER OUTPUT VARIABLE WHICH SPECIFIES THE C NUMBER OF GROUPS IN THE PARTITION OF THE COLUMNS OF A. C C IWA IS AN INTEGER WORK ARRAY OF LENGTH N. C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER IC,IP,IR,J,JCOL,JP C C INITIALIZATION BLOCK. C MAXGRP = 0 DO 10 JP = 1, N NGRP(JP) = N IWA(JP) = 0 10 CONTINUE C C BEGINNING OF ITERATION LOOP. C DO 60 J = 1, N JCOL = LIST(J) C C FIND ALL COLUMNS ADJACENT TO COLUMN JCOL. C C DETERMINE ALL POSITIONS (IR,JCOL) WHICH CORRESPOND C TO NON-ZEROES IN THE MATRIX. C DO 30 JP = JPNTR(JCOL), JPNTR(JCOL+1)-1 IR = INDROW(JP) C C FOR EACH ROW IR, DETERMINE ALL POSITIONS (IR,IC) C WHICH CORRESPOND TO NON-ZEROES IN THE MATRIX. C DO 20 IP = IPNTR(IR), IPNTR(IR+1)-1 IC = INDCOL(IP) C C ARRAY IWA MARKS THE GROUP NUMBERS OF THE C COLUMNS WHICH ARE ADJACENT TO COLUMN JCOL. C IWA(NGRP(IC)) = J 20 CONTINUE 30 CONTINUE C C ASSIGN THE SMALLEST UN-MARKED GROUP NUMBER TO JCOL. C DO 40 JP = 1, MAXGRP IF (IWA(JP) .NE. J) GO TO 50 40 CONTINUE MAXGRP = MAXGRP + 1 50 CONTINUE NGRP(JCOL) = JP 60 CONTINUE C C END OF ITERATION LOOP. C RETURN C C LAST CARD OF SUBROUTINE SEQ. C END SUBROUTINE SETR(M,N,INDROW,JPNTR,INDCOL,IPNTR,IWA) INTEGER M,N INTEGER INDROW(*),JPNTR(N+1),INDCOL(*),IPNTR(M+1),IWA(M) C ********** C C SUBROUTINE SETR C C GIVEN A COLUMN-ORIENTED DEFINITION OF THE SPARSITY PATTERN C OF AN M BY N MATRIX A, THIS SUBROUTINE DETERMINES A C ROW-ORIENTED DEFINITION OF THE SPARSITY PATTERN OF A. C C ON INPUT THE COLUMN-ORIENTED DEFINITION IS SPECIFIED BY C THE ARRAYS INDROW AND JPNTR. ON OUTPUT THE ROW-ORIENTED C DEFINITION IS SPECIFIED BY THE ARRAYS INDCOL AND IPNTR. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE SETR(M,N,INDROW,JPNTR,INDCOL,IPNTR,IWA) C C WHERE C C M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF ROWS OF A. C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF COLUMNS OF A. C C INDROW IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW C INDICES FOR THE NON-ZEROES IN THE MATRIX A. C C JPNTR IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW. C THE ROW INDICES FOR COLUMN J ARE C C INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1. C C NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO C ELEMENTS OF THE MATRIX A. C C INDCOL IS AN INTEGER OUTPUT ARRAY WHICH CONTAINS THE C COLUMN INDICES FOR THE NON-ZEROES IN THE MATRIX A. C C IPNTR IS AN INTEGER OUTPUT ARRAY OF LENGTH M + 1 WHICH C SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL. C THE COLUMN INDICES FOR ROW I ARE C C INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1. C C NOTE THAT IPNTR(1) IS SET TO 1 AND THAT IPNTR(M+1)-1 IS C THEN THE NUMBER OF NON-ZERO ELEMENTS OF THE MATRIX A. C C IWA IS AN INTEGER WORK ARRAY OF LENGTH M. C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER IR,JCOL,JP C C STORE IN ARRAY IWA THE COUNTS OF NON-ZEROES IN THE ROWS. C DO 10 IR = 1, M IWA(IR) = 0 10 CONTINUE DO 20 JP = 1, JPNTR(N+1)-1 IWA(INDROW(JP)) = IWA(INDROW(JP)) + 1 20 CONTINUE C C SET POINTERS TO THE START OF THE ROWS IN INDCOL. C IPNTR(1) = 1 DO 30 IR = 1, M IPNTR(IR+1) = IPNTR(IR) + IWA(IR) IWA(IR) = IPNTR(IR) 30 CONTINUE C C FILL INDCOL. C DO 50 JCOL = 1, N DO 40 JP = JPNTR(JCOL), JPNTR(JCOL+1)-1 IR = INDROW(JP) INDCOL(IWA(IR)) = JCOL IWA(IR) = IWA(IR) + 1 40 CONTINUE 50 CONTINUE RETURN C C LAST CARD OF SUBROUTINE SETR. C END SUBROUTINE SLO(N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,LIST, * MAXCLQ,IWA1,IWA2,IWA3,IWA4) INTEGER N,MAXCLQ INTEGER INDROW(*),JPNTR(N+1),INDCOL(*),IPNTR(*),NDEG(N), * LIST(N),IWA1(0:N-1),IWA2(N),IWA3(N),IWA4(N) C ********** C C SUBROUTINE SLO C C GIVEN THE SPARSITY PATTERN OF AN M BY N MATRIX A, THIS C SUBROUTINE DETERMINES THE SMALLEST-LAST ORDERING OF THE C COLUMNS OF A. C C THE SMALLEST-LAST ORDERING IS DEFINED FOR THE LOOPLESS C GRAPH G WITH VERTICES A(J), J = 1,2,...,N WHERE A(J) IS THE C J-TH COLUMN OF A AND WITH EDGE (A(I),A(J)) IF AND ONLY IF C COLUMNS I AND J HAVE A NON-ZERO IN THE SAME ROW POSITION. C C THE SMALLEST-LAST ORDERING IS DETERMINED RECURSIVELY BY C LETTING LIST(K), K = N,...,1 BE A COLUMN WITH LEAST DEGREE C IN THE SUBGRAPH SPANNED BY THE UN-ORDERED COLUMNS. C C NOTE THAT THE VALUE OF M IS NOT NEEDED BY SLO AND IS C THEREFORE NOT PRESENT IN THE SUBROUTINE STATEMENT. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE SLO(N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,LIST, C MAXCLQ,IWA1,IWA2,IWA3,IWA4) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF COLUMNS OF A. C C INDROW IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW C INDICES FOR THE NON-ZEROES IN THE MATRIX A. C C JPNTR IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW. C THE ROW INDICES FOR COLUMN J ARE C C INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1. C C NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO C ELEMENTS OF THE MATRIX A. C C INDCOL IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE C COLUMN INDICES FOR THE NON-ZEROES IN THE MATRIX A. C C IPNTR IS AN INTEGER INPUT ARRAY OF LENGTH M + 1 WHICH C SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL. C THE COLUMN INDICES FOR ROW I ARE C C INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1. C C NOTE THAT IPNTR(M+1)-1 IS THEN THE NUMBER OF NON-ZERO C ELEMENTS OF THE MATRIX A. C C NDEG IS AN INTEGER INPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE DEGREE SEQUENCE. THE DEGREE OF THE J-TH COLUMN C OF A IS NDEG(J). C C LIST IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE SMALLEST-LAST ORDERING OF THE COLUMNS OF A. THE J-TH C COLUMN IN THIS ORDER IS LIST(J). C C MAXCLQ IS AN INTEGER OUTPUT VARIABLE SET TO THE SIZE C OF THE LARGEST CLIQUE FOUND DURING THE ORDERING. C C IWA1,IWA2,IWA3, AND IWA4 ARE INTEGER WORK ARRAYS OF LENGTH N. C C SUBPROGRAMS CALLED C C FORTRAN-SUPPLIED ... MIN C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. AUGUST 1984. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER IC,IP,IR,JCOL,JP,MINDEG,NUMDEG,NUMORD C C INITIALIZATION BLOCK. C MINDEG = N DO 10 JP = 1, N IWA1(JP-1) = 0 IWA4(JP) = N LIST(JP) = NDEG(JP) MINDEG = MIN(MINDEG,NDEG(JP)) 10 CONTINUE C C CREATE A DOUBLY-LINKED LIST TO ACCESS THE DEGREES OF THE C COLUMNS. THE POINTERS FOR THE LINKED LIST ARE AS FOLLOWS. C C EACH UN-ORDERED COLUMN IC IS IN A LIST (THE DEGREE LIST) C OF COLUMNS WITH THE SAME DEGREE. C C IWA1(NUMDEG) IS THE FIRST COLUMN IN THE NUMDEG LIST C UNLESS IWA1(NUMDEG) = 0. IN THIS CASE THERE ARE C NO COLUMNS IN THE NUMDEG LIST. C C IWA2(IC) IS THE COLUMN BEFORE IC IN THE DEGREE LIST C UNLESS IWA2(IC) = 0. IN THIS CASE IC IS THE FIRST C COLUMN IN THIS DEGREE LIST. C C IWA3(IC) IS THE COLUMN AFTER IC IN THE DEGREE LIST C UNLESS IWA3(IC) = 0. IN THIS CASE IC IS THE LAST C COLUMN IN THIS DEGREE LIST. C C IF IC IS AN UN-ORDERED COLUMN, THEN LIST(IC) IS THE C DEGREE OF IC IN THE GRAPH INDUCED BY THE UN-ORDERED C COLUMNS. IF JCOL IS AN ORDERED COLUMN, THEN LIST(JCOL) C IS THE SMALLEST-LAST ORDER OF COLUMN JCOL. C DO 20 JP = 1, N NUMDEG = NDEG(JP) IWA2(JP) = 0 IWA3(JP) = IWA1(NUMDEG) IF (IWA1(NUMDEG) .GT. 0) IWA2(IWA1(NUMDEG)) = JP IWA1(NUMDEG) = JP 20 CONTINUE MAXCLQ = 0 NUMORD = N C C BEGINNING OF ITERATION LOOP. C 30 CONTINUE C C CHOOSE A COLUMN JCOL OF MINIMAL DEGREE MINDEG. C 40 CONTINUE JCOL = IWA1(MINDEG) IF (JCOL .GT. 0) GO TO 50 MINDEG = MINDEG + 1 GO TO 40 50 CONTINUE LIST(JCOL) = NUMORD C C MARK THE SIZE OF THE LARGEST CLIQUE C FOUND DURING THE ORDERING. C IF (MINDEG+1 .EQ. NUMORD .AND. MAXCLQ .EQ. 0) * MAXCLQ = NUMORD C C TERMINATION TEST. C NUMORD = NUMORD - 1 IF (NUMORD .EQ. 0) GO TO 80 C C DELETE COLUMN JCOL FROM THE MINDEG LIST. C IWA1(MINDEG) = IWA3(JCOL) IF (IWA3(JCOL) .GT. 0) IWA2(IWA3(JCOL)) = 0 C C FIND ALL COLUMNS ADJACENT TO COLUMN JCOL. C IWA4(JCOL) = 0 C C DETERMINE ALL POSITIONS (IR,JCOL) WHICH CORRESPOND C TO NON-ZEROES IN THE MATRIX. C DO 70 JP = JPNTR(JCOL), JPNTR(JCOL+1)-1 IR = INDROW(JP) C C FOR EACH ROW IR, DETERMINE ALL POSITIONS (IR,IC) C WHICH CORRESPOND TO NON-ZEROES IN THE MATRIX. C DO 60 IP = IPNTR(IR), IPNTR(IR+1)-1 IC = INDCOL(IP) C C ARRAY IWA4 MARKS COLUMNS WHICH ARE ADJACENT TO C COLUMN JCOL. C IF (IWA4(IC) .GT. NUMORD) THEN IWA4(IC) = NUMORD C C UPDATE THE POINTERS TO THE CURRENT DEGREE LISTS. C NUMDEG = LIST(IC) LIST(IC) = LIST(IC) - 1 MINDEG = MIN(MINDEG,LIST(IC)) C C DELETE COLUMN IC FROM THE NUMDEG LIST. C IF (IWA2(IC) .EQ. 0) THEN IWA1(NUMDEG) = IWA3(IC) ELSE IWA3(IWA2(IC)) = IWA3(IC) END IF IF (IWA3(IC) .GT. 0) IWA2(IWA3(IC)) = IWA2(IC) C C ADD COLUMN IC TO THE NUMDEG-1 LIST. C IWA2(IC) = 0 IWA3(IC) = IWA1(NUMDEG-1) IF (IWA1(NUMDEG-1) .GT. 0) IWA2(IWA1(NUMDEG-1)) = IC IWA1(NUMDEG-1) = IC END IF 60 CONTINUE 70 CONTINUE C C END OF ITERATION LOOP. C GO TO 30 80 CONTINUE C C INVERT THE ARRAY LIST. C DO 90 JCOL = 1, N IWA2(LIST(JCOL)) = JCOL 90 CONTINUE DO 100 JP = 1, N LIST(JP) = IWA2(JP) 100 CONTINUE RETURN C C LAST CARD OF SUBROUTINE SLO. C END SUBROUTINE SRTDAT(N,NNZ,INDROW,INDCOL,JPNTR,IWA) INTEGER N,NNZ INTEGER INDROW(NNZ),INDCOL(NNZ),JPNTR(N+1),IWA(N) C ********** C C SUBROUTINE SRTDAT C C GIVEN THE NON-ZERO ELEMENTS OF AN M BY N MATRIX A IN C ARBITRARY ORDER AS SPECIFIED BY THEIR ROW AND COLUMN C INDICES, THIS SUBROUTINE PERMUTES THESE ELEMENTS SO C THAT THEIR COLUMN INDICES ARE IN NON-DECREASING ORDER. C C ON INPUT IT IS ASSUMED THAT THE ELEMENTS ARE SPECIFIED IN C C INDROW(K),INDCOL(K), K = 1,...,NNZ. C C ON OUTPUT THE ELEMENTS ARE PERMUTED SO THAT INDCOL IS C IN NON-DECREASING ORDER. IN ADDITION, THE ARRAY JPNTR C IS SET SO THAT THE ROW INDICES FOR COLUMN J ARE C C INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1. C C NOTE THAT THE VALUE OF M IS NOT NEEDED BY SRTDAT AND IS C THEREFORE NOT PRESENT IN THE SUBROUTINE STATEMENT. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE SRTDAT(N,NNZ,INDROW,INDCOL,JPNTR,IWA) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF COLUMNS OF A. C C NNZ IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF NON-ZERO ELEMENTS OF A. C C INDROW IS AN INTEGER ARRAY OF LENGTH NNZ. ON INPUT INDROW C MUST CONTAIN THE ROW INDICES OF THE NON-ZERO ELEMENTS OF A. C ON OUTPUT INDROW IS PERMUTED SO THAT THE CORRESPONDING C COLUMN INDICES OF INDCOL ARE IN NON-DECREASING ORDER. C C INDCOL IS AN INTEGER ARRAY OF LENGTH NNZ. ON INPUT INDCOL C MUST CONTAIN THE COLUMN INDICES OF THE NON-ZERO ELEMENTS C OF A. ON OUTPUT INDCOL IS PERMUTED SO THAT THESE INDICES C ARE IN NON-DECREASING ORDER. C C JPNTR IS AN INTEGER OUTPUT ARRAY OF LENGTH N + 1 WHICH C SPECIFIES THE LOCATIONS OF THE ROW INDICES IN THE OUTPUT C INDROW. THE ROW INDICES FOR COLUMN J ARE C C INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1. C C NOTE THAT JPNTR(1) IS SET TO 1 AND THAT JPNTR(N+1)-1 C IS THEN NNZ. C C IWA IS AN INTEGER WORK ARRAY OF LENGTH N. C C SUBPROGRAMS CALLED C C FORTRAN-SUPPLIED ... MAX C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER I,J,K,L C C STORE IN ARRAY IWA THE COUNTS OF NON-ZEROES IN THE COLUMNS. C DO 10 J = 1, N IWA(J) = 0 10 CONTINUE DO 20 K = 1, NNZ IWA(INDCOL(K)) = IWA(INDCOL(K)) + 1 20 CONTINUE C C SET POINTERS TO THE START OF THE COLUMNS IN INDROW. C JPNTR(1) = 1 DO 30 J = 1, N JPNTR(J+1) = JPNTR(J) + IWA(J) IWA(J) = JPNTR(J) 30 CONTINUE K = 1 C C BEGIN IN-PLACE SORT. C 40 CONTINUE J = INDCOL(K) IF (K .GE. JPNTR(J)) THEN C C CURRENT ELEMENT IS IN POSITION. NOW EXAMINE THE C NEXT ELEMENT OR THE FIRST UN-SORTED ELEMENT IN C THE J-TH GROUP. C K = MAX(K+1,IWA(J)) ELSE C C CURRENT ELEMENT IS NOT IN POSITION. PLACE ELEMENT C IN POSITION AND MAKE THE DISPLACED ELEMENT THE C CURRENT ELEMENT. C L = IWA(J) IWA(J) = IWA(J) + 1 I = INDROW(K) INDROW(K) = INDROW(L) INDCOL(K) = INDCOL(L) INDROW(L) = I INDCOL(L) = J END IF IF (K .LE. NNZ) GO TO 40 RETURN C C LAST CARD OF SUBROUTINE SRTDAT. C END C ********** C C THIS IS A TEST PROGRAM FOR SUBROUTINES DSSM AND FDHS. C THE TEST DATA REPRESENTS SPECIAL TEST PROBLEMS. C C ********** INTEGER MAXNNZ,MAXN,LIWA PARAMETER (MAXNNZ = 38, MAXN = 12, LIWA = 72) C C LOGICAL OUTPUT UNIT IS ASSUMED TO BE NUMBER 6. C INTEGER NWRITE PARAMETER (NWRITE = 6) C INTEGER I,INFO,J,JP,MAXGRP,MAXROW,METHOD,MINGRP, * N,NNZ, NPAIRS,NUMGRP INTEGER INDCOL(MAXNNZ),INDROW(MAXNNZ),IPNTR(MAXN+1), * JPNTR(MAXN+1),LISTP(MAXN),NGRP(MAXN),IWA(LIWA) REAL ABSERR,DNSM,ENTRY1,ENTRY2,ERRIJ,FHEST,RELERR REAL ETA(MAXN),FHES(MAXNNZ),FHESD(MAXN),GVEC(MAXN), * X(MAXN),XD(MAXN) C C TEST FOR DSSM AND FDHS. C WRITE (NWRITE,1000) 10 CONTINUE C C DEFINITION OF SPARSITY PATTERN. C CALL TSTM(N,NPAIRS,INDROW,INDCOL,IWA) IF (N .EQ. 0) STOP DO 110 METHOD = 1, 2 C C CALL DSSM. C CALL DSSM(N,NPAIRS,INDROW,INDCOL,METHOD,LISTP,NGRP,MAXGRP, * MINGRP,INFO,IPNTR,JPNTR,IWA,LIWA) IF (INFO .LE. 0) THEN WRITE (NWRITE,2000) INFO STOP END IF C C STATISTICS FOR THE MATRIX. C MAXROW = 0 DO 20 I = 1, N MAXROW = MAX(MAXROW,IPNTR(I+1)-IPNTR(I)) 20 CONTINUE NNZ = IPNTR(N+1) - 1 DNSM = FLOAT(100*NNZ)/FLOAT(N*N) WRITE (NWRITE,3000) N,NNZ,DNSM,MAXROW,MINGRP,MAXGRP C C TEST FOR FDHS. C DO 30 J = 1, N X(J) = FLOAT(J)/FLOAT(N) 30 CONTINUE CALL FCN(N,X,INDROW,JPNTR,INDCOL,IPNTR,GVEC) C C APPROXIMATE THE HESSIAN MATRIX. C DO 60 NUMGRP = 1, MAXGRP DO 40 J = 1, N ETA(J) = 0.001 XD(J) = X(J) IF (NGRP(J) .EQ. NUMGRP) XD(J) = X(J) + ETA(J) 40 CONTINUE CALL FCN(N,XD,INDROW,JPNTR,INDCOL,IPNTR,FHESD) DO 50 I = 1, N FHESD(I) = FHESD(I) - GVEC(I) 50 CONTINUE CALL FDHS(N,INDROW,JPNTR,INDCOL,IPNTR, * LISTP,NGRP,MAXGRP,NUMGRP,ETA,FHESD,FHES,IWA) 60 CONTINUE C C TEST THE APPROXIMATION TO THE HESSIAN. C ABSERR = 0.0 RELERR = 0.0 ENTRY1 = 0.0 ENTRY2 = 0.0 DO 80 J = 1, N DO 70 JP = JPNTR(J), JPNTR(J+1)-1 I = INDROW(JP) FHEST = I + J FHEST = 1/FHEST IF (I .EQ. J) FHEST = 2*FHEST ERRIJ = FHES(JP) - FHEST IF (ABS(ERRIJ) .GT. ABSERR) THEN ENTRY1 = FHEST ABSERR = ABS(ERRIJ) END IF IF (ABS(ERRIJ) .GT. RELERR*ABS(FHEST)) THEN ENTRY2 = FHEST IF (FHEST .NE. 0.0) * RELERR = ABS(ERRIJ)/ABS(FHEST) END IF 70 CONTINUE 80 CONTINUE WRITE (NWRITE,4000) ABSERR,ENTRY1,RELERR,ENTRY2 DO 100 J = 1, N DO 90 JP = JPNTR(J), JPNTR(J+1)-1 INDCOL(JP) = J 90 CONTINUE 100 CONTINUE 110 CONTINUE GO TO 10 C C FORMAT STATEMENTS. C 1000 FORMAT (// ' TESTS FOR DSSM - TEST PROBLEMS' // * ' STATISTICS GENERATED '// * ' N - NUMBER OF COLUMNS '/ * ' NNZ - NUMBER OF NONZERO ELEMENTS'/ * ' DNSM - MATRIX DENSITY (PERCENTAGE)'/ * ' MAXROW - MAXIMUM NUMBER OF NONZEROES IN ANY ROW'/ * ' MINGRP - LOWER BOUND ON NUMBER OF GROUPS'/ * ' MAXGRP - NUMBER OF GROUPS DETERMINED BY DSSM'///) 2000 FORMAT (// ' *** MISTAKE IN INPUT DATA, INFO IS ***',I6) 3000 FORMAT (// 3X,'N',6X,'NNZ',5X,'DNSM',5X, * 'MAXROW',4X,'MINGRP',4X,'MAXGRP'// * 2(I5,3X),F6.2,4X,3(I5,5X)) 4000 FORMAT (// ' LARGEST ABSOLUTE ERROR ',E10.2,5X,'ENTRY',E10.2, * / ' LARGEST RELATIVE ERROR ',E10.2,5X,'ENTRY',E10.2) C C LAST CARD OF MAIN PROGRAM. C END SUBROUTINE FCN(N,X,INDROW,JPNTR,INDCOL,IPNTR,GVEC) INTEGER N INTEGER INDROW(*),JPNTR(N+1),INDCOL(*),IPNTR(N+1) REAL X(N),GVEC(N) C C GRADIENT SUBROUTINE FOR TESTING FDHS. C INTEGER I,IP REAL SUM C DO 30 I = 1, N SUM = 0.0 DO 10 IP = IPNTR(I), IPNTR(I+1)-1 SUM = SUM + X(INDCOL(IP))/(I+INDCOL(IP)) 10 CONTINUE DO 20 JP = JPNTR(I), JPNTR(I+1)-1 SUM = SUM + X(INDROW(JP))/(I+INDROW(JP)) 20 CONTINUE GVEC(I) = SUM 30 CONTINUE RETURN C C LAST CARD OF SUBROUTINE FCN. C END SUBROUTINE TSTM(N,NNZ,INDROW,INDCOL,IWA) INTEGER N,NNZ INTEGER INDROW(*),INDCOL(*),IWA(*) C ********** C C SUBROUTINE TSTMAT C C ********** INTEGER I,J,MAXDEG,NCOL C C LOGICAL INPUT UNIT IS ASSUMED TO BE NUMBER 5. C INTEGER NREAD PARAMETER ( = 5) C READ (NREAD,40) N,MAXDEG NNZ = 0 DO 30 J = 1, N READ (NREAD,40) NCOL,(IWA(I),I=1,MAXDEG) DO 10 I = 1, MAXDEG IF (IWA(I) .EQ. 0) GO TO 20 IF (IWA(I) .GE. NCOL) THEN NNZ = NNZ + 1 INDROW(NNZ) = IWA(I) INDCOL(NNZ) = NCOL END IF 10 CONTINUE 20 CONTINUE 30 CONTINUE 40 FORMAT (16I5) RETURN C C LAST CARD OF SUBROUTINE TSTMAT. C END 1 1 1 1 4 1 1 1 2 2 3 3 4 4 7 3 1 1 2 0 2 1 2 3 3 2 3 4 4 3 4 5 5 4 5 6 6 5 6 7 7 6 7 0 8 5 1 1 2 3 0 2 1 2 3 4 0 3 1 2 3 4 5 4 2 3 4 5 6 5 3 4 5 6 7 6 4 5 6 7 8 7 5 6 7 8 0 8 6 7 8 0 11 7 1 1 2 3 4 0 2 1 2 3 4 5 0 3 1 2 3 4 5 6 0 4 1 2 3 4 5 6 7 5 2 3 4 5 6 7 8 6 3 4 5 6 7 8 9 7 4 5 6 7 8 9 10 8 5 6 7 8 9 10 11 9 6 7 8 9 10 11 0 10 7 8 9 10 11 0 11 8 9 10 11 0 12 5 1 1 2 3 0 2 1 2 4 5 0 3 1 3 6 7 8 4 2 4 9 10 11 5 2 5 0 6 3 6 0 7 3 7 0 8 3 8 0 9 4 9 0 10 4 10 0 11 4 11 12 0 12 11 12 0 9 7 1 1 4 5 6 7 8 9 2 2 4 5 6 7 8 9 3 3 4 5 6 7 8 9 4 4 1 2 3 0 5 5 1 2 3 0 6 6 1 2 3 0 7 7 1 2 3 0 8 8 1 2 3 0 9 9 1 2 3 0 5 5 1 1 2 3 4 5 2 1 2 3 4 5 3 1 2 3 4 5 4 1 2 3 4 5 5 1 2 3 4 5 6 4 1 1 2 3 4 2 1 2 3 5 3 1 2 3 6 4 1 4 0 5 2 5 0 6 3 6 0 5 5 1 1 5 0 2 2 5 0 3 3 5 0 4 4 5 0 5 1 2 3 4 5 0 0 0 0 0 7 4 1 1 5 6 7 2 2 3 4 0 3 2 3 4 0 4 2 3 4 0 5 1 5 0 6 1 6 0 7 1 7 0 0 0 0 0 0