C ALGORITHM 618, COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.10, NO. 3, C SEP., 1984, P. 346-347. PROGRAM TEST C ********** C C THIS IS A TEST PROGRAM FOR SUBROUTINES DSM AND FDJS. C THE TEST DATA REPRESENTS A NEUTRON KINETICS PROBLEM. C C ********** INTEGER I,INFO,IP,J,JP,L,LIWA,M, * MAXGRP,MAXROW,MINGRP,MINROW,N,NNZ,NUMGRP,NWRITE INTEGER INDCOL(6000),INDROW(6000),IPNTR(1201),JPNTR(1201), * NGRP(1200),IWA(7200) LOGICAL COL REAL DNSM,ERRIJ,ERRMAX,FJACT,FJACTR,SUM REAL D(1200),FJAC(6000),FJACD(1200),FVEC(1200),X(1200),XD(1200) C C LOGICAL OUTPUT UNIT IS ASSUMED TO BE NUMBER 6. C DATA NWRITE /6/ C LIWA = 7200 COL = .TRUE. C C TEST FOR DSM AND FDJS. C WRITE (NWRITE,1000) DO 130 N = 300, 1200, 300 WRITE (NWRITE,2000) C C DEFINITION OF SPARSITY PATTERN. C M = N L = N/3 NNZ = 0 DO 10 J = 1, N NNZ = NNZ + 1 INDROW(NNZ) = J INDCOL(NNZ) = J IF (MOD(J,L) .NE. 0) THEN NNZ = NNZ + 1 INDROW(NNZ) = J + 1 INDCOL(NNZ) = J END IF IF (J .LE. 2*L) THEN NNZ = NNZ + 1 INDROW(NNZ) = J + L INDCOL(NNZ) = J IF (MOD(J,L) .NE. 1) THEN NNZ = NNZ + 1 INDROW(NNZ) = J - 1 INDCOL(NNZ) = J END IF END IF NNZ = NNZ + 1 IF (J .GT. L) THEN INDROW(NNZ) = J - L ELSE INDROW(NNZ) = J + 2*L END IF INDCOL(NNZ) = J 10 CONTINUE C C CALL DSM. C CALL DSM(M,N,NNZ,INDROW,INDCOL,NGRP,MAXGRP,MINGRP, * INFO,IPNTR,JPNTR,IWA,LIWA) IF (INFO .LE. 0) WRITE (NWRITE,4000) INFO C C STATISTICS FOR THE MATRIX. C MAXROW = 0 MINROW = N DO 20 I = 1, M MAXROW = MAX(MAXROW,IPNTR(I+1)-IPNTR(I)) MINROW = MIN(MINROW,IPNTR(I+1)-IPNTR(I)) 20 CONTINUE DNSM = FLOAT(100*NNZ)/FLOAT(M*N) WRITE (NWRITE,3000) N,NNZ,DNSM,MINROW,MAXROW,MINGRP,MAXGRP C C TEST FOR FDJS. C DO 30 J = 1, N X(J) = FLOAT(J)/FLOAT(N) 30 CONTINUE CALL FCN(N,X,INDCOL,IPNTR,FVEC) C C APPROXIMATE THE JACOBIAN MATRIX. C DO 60 NUMGRP = 1, MAXGRP DO 40 J = 1, N D(J) = 0.0 IF (NGRP(J) .EQ. NUMGRP) D(J) = 0.001 XD(J) = X(J) + D(J) 40 CONTINUE CALL FCN(N,XD,INDCOL,IPNTR,FJACD) DO 50 I = 1, M FJACD(I) = FJACD(I) - FVEC(I) 50 CONTINUE IF (COL) THEN CALL FDJS(M,N,COL,INDROW,JPNTR,NGRP,NUMGRP,D,FJACD,FJAC) ELSE CALL FDJS(M,N,COL,INDCOL,IPNTR,NGRP,NUMGRP,D,FJACD,FJAC) END IF 60 CONTINUE C C TEST THE APPROXIMATION TO THE JACOBIAN. C ERRMAX = 0.0 IF (COL) THEN C C TEST FOR THE COLUMN-ORIENTED DEFINITION OF C THE SPARSITY PATTERN. C DO 90 J = 1, N DO 80 JP = JPNTR(J), JPNTR(J+1)-1 I = INDROW(JP) SUM = 0.0 DO 70 IP = IPNTR(I), IPNTR(I+1)-1 SUM = SUM + X(INDCOL(IP)) 70 CONTINUE SUM = SUM + X(I) FJACT = 1.0 + 2.0*SUM IF (I .EQ. J) FJACT = 2.0*FJACT ERRIJ = FJAC(JP) - FJACT IF (FJACT .NE. 0.0) ERRIJ = ERRIJ/FJACT ERRMAX = MAX(ERRMAX,ABS(ERRIJ)) 80 CONTINUE 90 CONTINUE ELSE C C TEST FOR THE ROW-ORIENTED DEFINITION OF C THE SPARSITY PATTERN. C DO 120 I = 1, M SUM = 0.0 DO 100 IP = IPNTR(I), IPNTR(I+1)-1 SUM = SUM + X(INDCOL(IP)) 100 CONTINUE SUM = SUM + X(I) FJACTR = 1.0 + 2.0*SUM DO 110 IP = IPNTR(I), IPNTR(I+1)-1 J = INDCOL(IP) FJACT = FJACTR IF (I .EQ. J) FJACT = 2.0*FJACT ERRIJ = FJAC(IP) - FJACT IF (FJACT .NE. 0.0) ERRIJ = ERRIJ/FJACT ERRMAX = MAX(ERRMAX,ABS(ERRIJ)) 110 CONTINUE 120 CONTINUE END IF WRITE (NWRITE,5000) ERRMAX COL = .NOT. COL 130 CONTINUE STOP C C FORMAT STATEMENTS. C 1000 FORMAT(// ' TESTS FOR DSM AND FDJS - NEUTRON KINETICS PROBLEM'// * ' STATISTICS GENERATED '// * ' N - NUMBER OF COLUMNS '/ * ' NNZ - NUMBER OF NON-ZERO ELEMENTS'/ * ' DNSM - MATRIX DENSITY (PERCENTAGE)'/ * ' MINROW - MINIMUM NUMBER OF NON-ZEROS IN ANY ROW'/ * ' MAXROW - MAXIMUM NUMBER OF NON-ZEROS IN ANY ROW'/ * ' MINGRP - LOWER BOUND ON NUMBER OF GROUPS'/ * ' MAXGRP - NUMBER OF GROUPS DETERMINED BY DSM'//) 2000 FORMAT(// 3X,'N',6X,'NNZ',5X,'DNSM',5X, * 'MINROW',4X,'MAXROW',4X,'MINGRP',4X,'MAXGRP'//) 3000 FORMAT(2(I5,3X),F6.2,4X,4(I5,5X)) 4000 FORMAT(// ' *** MISTAKE IN INPUT DATA, INFO IS ***',I6) 5000 FORMAT(// ' LARGEST RELATIVE ERROR OF APPROXIMATION IS',E10.2) END SUBROUTINE FCN(N,X,INDCOL,IPNTR,FVEC) INTEGER N INTEGER INDCOL(*),IPNTR(N+1) REAL X(N),FVEC(N) C C FUNCTION SUBROUTINE FOR TESTING FDJS. C INTEGER I,IP REAL SUM DO 20 I = 1, N SUM = 0.0 DO 10 IP = IPNTR(I), IPNTR(I+1)-1 SUM = SUM + X(INDCOL(IP)) 10 CONTINUE SUM = SUM + X(I) FVEC(I) = SUM*(1.0 + SUM) + 1.0 20 CONTINUE RETURN C C LAST CARD OF SUBROUTINE FCN. C END SUBROUTINE DSM(M,N,NPAIRS,INDROW,INDCOL,NGRP,MAXGRP,MINGRP, * INFO,IPNTR,JPNTR,IWA,LIWA) INTEGER M,N,NPAIRS,MAXGRP,MINGRP,INFO,LIWA INTEGER INDROW(NPAIRS),INDCOL(NPAIRS),NGRP(N), * IPNTR(M+1),JPNTR(N+1),IWA(LIWA) C ********** C C SUBROUTINE DSM C C THE PURPOSE OF DSM IS TO DETERMINE AN OPTIMAL OR NEAR- C OPTIMAL CONSISTENT PARTITION OF THE COLUMNS OF A SPARSE C M BY N MATRIX A. C C THE SPARSITY PATTERN OF THE MATRIX A IS SPECIFIED BY C THE ARRAYS INDROW AND INDCOL. ON INPUT THE INDICES C FOR THE NON-ZERO ELEMENTS OF A ARE C C INDROW(K),INDCOL(K), K = 1,2,...,NPAIRS. C C THE (INDROW,INDCOL) PAIRS MAY BE SPECIFIED IN ANY ORDER. C DUPLICATE INPUT PAIRS ARE PERMITTED, BUT THE SUBROUTINE C ELIMINATES THEM. C C THE SUBROUTINE PARTITIONS THE COLUMNS OF A INTO GROUPS C SUCH THAT COLUMNS IN THE SAME GROUP DO NOT HAVE A C NON-ZERO IN THE SAME ROW POSITION. A PARTITION OF THE C COLUMNS OF A WITH THIS PROPERTY IS CONSISTENT WITH THE C DIRECT DETERMINATION OF A. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE DSM(M,N,NPAIRS,INDROW,INDCOL,NGRP,MAXGRP,MINGRP, C INFO,IPNTR,JPNTR,IWA,LIWA) 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 NPAIRS IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE C NUMBER OF (INDROW,INDCOL) PAIRS USED TO DESCRIBE THE C SPARSITY 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 OF A. C ON OUTPUT INDROW IS PERMUTED SO THAT THE CORRESPONDING C COLUMN INDICES ARE IN NON-DECREASING ORDER. THE COLUMN C INDICES CAN BE RECOVERED 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 OF C A. ON OUTPUT INDCOL IS PERMUTED SO THAT THE CORRESPONDING C ROW INDICES ARE IN NON-DECREASING ORDER. THE ROW INDICES C CAN BE RECOVERED FROM THE ARRAY IPNTR. 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 MINGRP IS AN INTEGER OUTPUT VARIABLE WHICH SPECIFIES A LOWER C BOUND FOR THE NUMBER OF GROUPS IN ANY CONSISTENT PARTITION C OF THE COLUMNS OF A. C C INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS. FOR C NORMAL TERMINATION INFO = 1. IF M, N, OR NPAIRS IS NOT C POSITIVE OR LIWA IS LESS THAN MAX(M,6*N), THEN INFO = 0. C IF THE K-TH ELEMENT OF INDROW IS NOT AN INTEGER BETWEEN C 1 AND M OR THE K-TH ELEMENT OF INDCOL IS NOT AN INTEGER C BETWEEN 1 AND N, THEN INFO = -K. 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(M+1)-1 IS THEN THE NUMBER OF NON-ZERO C ELEMENTS 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 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 C MAX(M,6*N). C C SUBPROGRAMS CALLED C C MINPACK-SUPPLIED ... DEGR,IDO,NUMSRT,SEQ,SETR,SLO,SRTDAT 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,IR,J,JP,K,MAXCLQ,NNZ,NUMGRP C C CHECK THE INPUT DATA. C INFO = 0 IF (M .LT. 1 .OR. N .LT. 1 .OR. NPAIRS .LT. 1 .OR. * LIWA .LT. MAX(M,6*N)) RETURN DO 10 K = 1, NPAIRS INFO = -K IF (INDROW(K) .LT. 1 .OR. INDROW(K) .GT. M .OR. * INDCOL(K) .LT. 1 .OR. INDCOL(K) .GT. N) RETURN 10 CONTINUE INFO = 1 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 C NON-ZERO ELEMENTS OF A. C DO 20 I = 1, M IWA(I) = 0 20 CONTINUE NNZ = 1 DO 40 J = 1, N K = NNZ DO 30 JP = JPNTR(J), JPNTR(J+1)-1 IR = INDROW(JP) IF (IWA(IR) .NE. J) THEN INDROW(NNZ) = IR NNZ = NNZ + 1 IWA(IR) = J END IF 30 CONTINUE JPNTR(J) = K 40 CONTINUE JPNTR(N+1) = NNZ C C EXTEND THE DATA STRUCTURE TO ROWS. C CALL SETR(M,N,INDROW,JPNTR,INDCOL,IPNTR,IWA) C C DETERMINE A LOWER BOUND FOR THE NUMBER OF GROUPS. C MINGRP = 0 DO 50 I = 1, M MINGRP = MAX(MINGRP,IPNTR(I+1)-IPNTR(I)) 50 CONTINUE C C DETERMINE THE DEGREE SEQUENCE FOR THE INTERSECTION C GRAPH OF THE COLUMNS OF A. 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 A 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),NGRP,MAXGRP, * IWA(N+1)) MINGRP = MAX(MINGRP,MAXCLQ) C C EXIT IF THE SMALLEST-LAST ORDERING IS OPTIMAL. C IF (MAXGRP .EQ. MINGRP) RETURN C C COLOR THE INTERSECTION GRAPH OF THE COLUMNS OF A C WITH THE INCIDENCE-DEGREE (ID) ORDERING. C CALL IDO(M,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)) MINGRP = MAX(MINGRP,MAXCLQ) C C RETAIN THE BETTER OF THE TWO ORDERINGS SO FAR. C IF (NUMGRP .LT. MAXGRP) THEN MAXGRP = NUMGRP DO 60 J = 1, N NGRP(J) = IWA(J) 60 CONTINUE C C EXIT IF THE INCIDENCE-DEGREE ORDERING IS OPTIMAL. C IF (MAXGRP .EQ. MINGRP) RETURN END IF C C COLOR THE INTERSECTION GRAPH OF THE COLUMNS OF A C WITH THE LARGEST-FIRST (LF) ORDERING. C CALL NUMSRT(N,N-1,IWA(5*N+1),-1,IWA(4*N+1),IWA(2*N+1),IWA(N+1)) CALL SEQ(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(4*N+1),IWA(1),NUMGRP, * IWA(N+1)) C C RETAIN THE BEST OF THE THREE ORDERINGS AND EXIT. C IF (NUMGRP .LT. MAXGRP) THEN MAXGRP = NUMGRP DO 70 J = 1, N NGRP(J) = IWA(J) 70 CONTINUE END IF RETURN C C LAST CARD OF SUBROUTINE DSM. 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. JULY 1983. 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 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 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 NUMORD = NUMORD + 1 C C TERMINATION TEST. C 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. JULY 1983. 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 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 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 NUMORD = NUMORD - 1 C C TERMINATION TEST. C 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 SUBROUTINE FDJS(M,N,COL,IND,NPNTR,NGRP,NUMGRP,D,FJACD,FJAC) INTEGER M,N,NUMGRP INTEGER IND(*),NPNTR(*),NGRP(N) REAL D(N),FJACD(M),FJAC(*) LOGICAL COL C ********** C C SUBROUTINE FDJS C C GIVEN A CONSISTENT PARTITION OF THE COLUMNS OF AN M BY N C JACOBIAN MATRIX INTO GROUPS, THIS SUBROUTINE COMPUTES C APPROXIMATIONS TO THOSE COLUMNS IN A GIVEN GROUP. THE C APPROXIMATIONS ARE STORED INTO EITHER A COLUMN-ORIENTED C OR A ROW-ORIENTED PATTERN. C C A PARTITION IS CONSISTENT IF THE COLUMNS IN ANY GROUP C DO NOT HAVE A NON-ZERO IN THE SAME ROW POSITION. C C APPROXIMATIONS TO THE COLUMNS OF THE JACOBIAN MATRIX IN A C GIVEN GROUP CAN BE OBTAINED BY SPECIFYING A DIFFERENCE C PARAMETER ARRAY D WITH D(JCOL) NON-ZERO IF AND ONLY IF C JCOL IS A COLUMN IN THE GROUP, AND AN APPROXIMATION TO C JAC*D WHERE JAC DENOTES THE JACOBIAN MATRIX OF A MAPPING F. C C D CAN BE DEFINED WITH THE FOLLOWING SEGMENT OF CODE. C C DO 10 JCOL = 1, N C D(JCOL) = 0.0 C IF (NGRP(JCOL) .EQ. NUMGRP) D(JCOL) = ETA(JCOL) C 10 CONTINUE C C IN THE ABOVE CODE NUMGRP IS THE GIVEN GROUP NUMBER, C NGRP(JCOL) IS THE GROUP NUMBER OF COLUMN JCOL, AND C ETA(JCOL) IS THE DIFFERENCE PARAMETER USED TO C APPROXIMATE COLUMN JCOL OF THE JACOBIAN MATRIX. C SUITABLE VALUES FOR THE ARRAY ETA MUST BE PROVIDED. C C AS MENTIONED ABOVE, AN APPROXIMATION TO JAC*D MUST C ALSO BE PROVIDED. FOR EXAMPLE, THE APPROXIMATION C C F(X+D) - F(X) C C CORRESPONDS TO THE FORWARD DIFFERENCE FORMULA AT X. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE FDJS(M,N,COL,IND,NPNTR,NGRP,NUMGRP,D,FJACD,FJAC) C C WHERE C C M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF ROWS OF THE JACOBIAN MATRIX. C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF COLUMNS OF THE JACOBIAN MATRIX. C C COL IS A LOGICAL INPUT VARIABLE. IF COL IS SET TRUE, THEN THE C JACOBIAN APPROXIMATIONS ARE STORED INTO A COLUMN-ORIENTED C PATTERN. IF COL IS SET FALSE, THEN THE JACOBIAN C APPROXIMATIONS ARE STORED INTO A ROW-ORIENTED PATTERN. C C IND IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW C INDICES FOR THE NON-ZEROES IN THE JACOBIAN MATRIX C IF COL IS TRUE, AND CONTAINS THE COLUMN INDICES FOR C THE NON-ZEROES IN THE JACOBIAN MATRIX IF COL IS FALSE. C C NPNTR IS AN INTEGER INPUT ARRAY WHICH SPECIFIES THE C LOCATIONS OF THE ROW INDICES IN IND IF COL IS TRUE, AND C SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN IND IF C COL IS FALSE. IF COL IS TRUE, THE INDICES FOR COLUMN J ARE C C IND(K), K = NPNTR(J),...,NPNTR(J+1)-1. C C IF COL IS FALSE, THE INDICES FOR ROW I ARE C C IND(K), K = NPNTR(I),...,NPNTR(I+1)-1. C C NOTE THAT NPNTR(N+1)-1 IF COL IS TRUE, OR NPNTR(M+1)-1 C IF COL IS FALSE, IS THEN THE NUMBER OF NON-ZERO ELEMENTS C OF THE JACOBIAN MATRIX. C C NGRP IS AN INTEGER INPUT ARRAY OF LENGTH N WHICH SPECIFIES C THE PARTITION OF THE COLUMNS OF THE JACOBIAN MATRIX. C COLUMN JCOL BELONGS TO GROUP NGRP(JCOL). C C NUMGRP IS A POSITIVE INTEGER INPUT VARIABLE SET TO A GROUP C NUMBER IN THE PARTITION. THE COLUMNS OF THE JACOBIAN C MATRIX IN THIS GROUP ARE TO BE ESTIMATED ON THIS CALL. C C D IS AN INPUT ARRAY OF LENGTH N WHICH CONTAINS THE C DIFFERENCE PARAMETER VECTOR FOR THE ESTIMATE OF C THE JACOBIAN MATRIX COLUMNS IN GROUP NUMGRP. C C FJACD IS AN INPUT ARRAY OF LENGTH M WHICH CONTAINS C AN APPROXIMATION TO THE DIFFERENCE VECTOR JAC*D, C WHERE JAC DENOTES THE JACOBIAN MATRIX. C C FJAC IS AN OUTPUT ARRAY OF LENGTH NNZ, WHERE NNZ IS THE C NUMBER OF ITS NON-ZERO ELEMENTS. AT EACH CALL OF FDJS, C FJAC IS UPDATED TO INCLUDE THE NON-ZERO ELEMENTS OF THE C JACOBIAN MATRIX FOR THOSE COLUMNS IN GROUP NUMGRP. FJAC C SHOULD NOT BE ALTERED BETWEEN SUCCESSIVE CALLS TO FDJS. C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983. C THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE' C C ********** INTEGER IP,IROW,JCOL,JP C C COMPUTE ESTIMATES OF JACOBIAN MATRIX COLUMNS IN GROUP C NUMGRP. THE ARRAY FJACD MUST CONTAIN AN APPROXIMATION C TO JAC*D, WHERE JAC DENOTES THE JACOBIAN MATRIX AND D C IS A DIFFERENCE PARAMETER VECTOR WITH D(JCOL) NON-ZERO C IF AND ONLY IF JCOL IS A COLUMN IN GROUP NUMGRP. C IF (COL) THEN C C COLUMN ORIENTATION. C DO 20 JCOL = 1, N IF (NGRP(JCOL) .EQ. NUMGRP) THEN DO 10 JP = NPNTR(JCOL), NPNTR(JCOL+1)-1 IROW = IND(JP) FJAC(JP) = FJACD(IROW)/D(JCOL) 10 CONTINUE END IF 20 CONTINUE ELSE C C ROW ORIENTATION. C DO 50 IROW = 1, M DO 30 IP = NPNTR(IROW), NPNTR(IROW+1)-1 JCOL = IND(IP) IF (NGRP(JCOL) .EQ. NUMGRP) THEN FJAC(IP) = FJACD(IROW)/D(JCOL) GO TO 40 END IF 30 CONTINUE 40 CONTINUE 50 CONTINUE END IF RETURN C C LAST CARD OF SUBROUTINE FDJS. C END