C COPYRIGHT (c) 1997 Council for the Central Laboratory * of the Research Councils C######DATE 3 September 1997 C April 2001: call to MC49 changed to MC59 to make routine threadsafe SUBROUTINE MC54A(ICNTL,TITLE,KEY,M,N,NE,IP,IND,VALUE,IW,INFO) INTEGER ICNTL(10) CHARACTER TITLE*72,KEY*8 INTEGER M,N,NE,IP(*),IND(*) REAL VALUE(*) INTEGER IW(*),INFO(5) C ================================================================ C Code for writing a sparse matrix in Rutherford-Boeing format, C input matrix can be assembled (row ptr/col index or coordinate C form) or unassembled (finite-element form). C For matrices with complex entries MC54AC/AZ should be used. C For matrices with integer entries MC54AI should be used. C ================================================================ C There follows a quick resume of the arguments ... full details follow. C ICNTL(1) = output unit C ICNTL(2) = precision of real data C ICNTL(3) = numerical values or only pattern C ICNTL(4) = shape/symmetry C ICNTL(5) = type of input C TITLE = descriptive title for matrix C KEY = matrix identifier C M = number of rows or maximum index C N = number of columns or number of elements C NE = number of nonzeros or elemental indices C IP = pointer array C IND = index array C VALUE = numerical values C IW = work array C INFO = information and error flags C C ICNTL is an INTEGER array of length 10 that must be set by the user. C It holds information on the input matrix and output file and controls C the action of the subroutine. This argument is not altered by the C subroutine. C ICNTL(1) must be set by the user to the stream number for the output C file on which the Rutherford-Boeing format will be written. C ICNTL(2) must be set by the user to indicate the precision of the C output reals. It should be set to the number of significant C decimal digits in the numerical value of the entries. Values C of ICNTL(2) less than or equal to one or greater than 17 C are treated as if they were 17. ICNTL(2) is not accessed C if ICNTL(3)=0. C ICNTL(3) must be set by the user to indicate whether only the C pattern is being supplied (ICNTL(3)=0) or whether real C values are supplied (ICNTL(3) not equal to 0). C ICNTL(4) must be set by the user to indicate whether matrix is C symmetric (ICNTL(4)=0), skew symmetric (ICNTL(4)=1), C rectangular (ICNTL(4)=2), or unsymmetric (ICNTL(4)=3). C A value of ICNTL(4)=1 is invalid if ICNTL(3)=0. C An error flag is set if ICNTL(4) < 0 or > 3. C ICNTL(5) must be set by the user to indicate the format of the input C matrix: ICNTL(5)=0 for an assembled matrix in C column-oriented form, ICNTL(5)=1 for an assembled matrix in C coordinate form, and ICNTL(5)=2 for an unassembled C finite-element matrix. C An error flag is set if ICNTL(5) < 0 or > 2. C ICNTL(6) to ICNTL(10) are not accessed by the routine. C C TITLE is a CHARACTER*72 variable that must be set by the user to the C text for the title line for the matrix. This argument is not altered C by the subroutine. C C KEY is a CHARACTER*8 variable that must be set by the user to the C identifying key for the matrix. If the keepers of the Collection C find that this is identical to an already existing matrix, it will C be changed after consultation with submitter. C This argument is not altered by the subroutine. C C M is an INTEGER variable that must be set by the user. If the matrix C is assembled, M is the number of rows of the matrix; if the matrix is C an unassembled finite-element matrix, M is the value of the largest C index used. C This argument is not altered by the subroutine. C C N is an INTEGER variable that must be set by the user. If the matrix C is assembled, N is the number of columns of the matrix; if the matrix C is an unassembled finite-element matrix, N is the number of elements. C This argument is not altered by the subroutine. C C NE is an INTEGER variable that must be set by the user. If the matrix C is assembled, NE is the number of entries in the matrix; if the matrix C is an unassembled finite-element matrix, NE is the number of indices C in the element lists. C This argument is not altered by the subroutine. C C IP is an INTEGER array of size MAX(M,N)+1 or 2*N+1 if input is a C rectangular elemental matrix. If the matrix is stored in C column-oriented form, IP must be set by the user so that IP(J), C J=1, ..., N holds the position in IND of the first entry in column J C and IP(N+1)-1 holds the position of the last entry in column N. C For unassembled finite-element problems (ICNTL(5) not equal to 0 C or 1), IP(I) must point to the start of the list of indices for C element I in array IND. C If elements are rectangular, then the row indices for element I C are in positions IP(2*I-1) to IP(2*I)-1 of IND and the column C indices for element I are in positions IP(2*I) to IP(2*I+1)-1. C In these cases, this argument is not altered by the subroutine. C If the matrix is held in coordinate form, IP need not be set on entry C and, on exit, IP(J) will be set to the position in IND of the first C row index in column J in the reordered matrix. C C IND is an INTEGER array whose first NE positions must be set by C the user. If the matrix is assembled, IND must hold the row indices C of the entries. If storage is column-oriented, these must be ordered C by columns. If storage is in coordinate form, the corresponding column C indices must be in positions IND(NE+I), I=1, ... NE. In this case, C the matrix entries can be in any order. For a symmetric matrix, only C the lower triangle should be stored. If the matrix is an unassembled C finite-element matrix, IND must hold the list of indices of the C variables in each element in turn. IND is not altered by the C subroutine, unless input is in coordinate form (ICNTL(5)=1), in which C case it will hold row indices of the reordered matrix. C C VALUE is a REAL array. If ICNTL(3) is equal to 1, it C must be set by the user to contain the numerical values of the entries. C For an assembled matrix, the first NE positions must contain the C values of the entries whose row indices are in corresponding positions C in IND. For an unassembled finite-element matrix, the first entries C of VALUE must contain contiguously the numerical values of each C element. If these are symmetric only the lower triangle should be C held ordered by columns. VALUE is not altered by the subroutine, C unless matrix is input in coordinate form (ICNTL(5)=1), in which C case it will hold the values for the reordered matrix. C IW is an INTEGER work array of length MAX(M,N)+1 that is only accessed C if the matrix is input in coordinate form (ICNTL(5)=1). C C INFO is an INTEGER array that need not be set by the user. On exit, C INFO(1) will be zero if no error is detected, and negative otherwise. C Possible nonzero values are: C -1 M <= 0. INFO(2) holds value of M. C -2 N <= 0. INFO(2) holds value of N. C -3 NE <= 0. INFO(2) holds value of NE. C -4 ICNTL(4) out of range. INFO(2) holds value of ICNTL(4). C -5 ICNTL(5) out of range. INFO(2) holds value of ICNTL(5). C -6 ICNTL(4) = 1 but ICNTL(3) = 0. INFO(2) = 0. C -7 Matrix was submitted in column-oriented form with IP(I+1) < IP(I) C for some I, 1<= I <= N. INFO(2) holds the value of I. C -8 Index out-of-range. INFO(2) holds index of component. C -9 Matrix was submitted in column-oriented form with rows not in order C within each column. INFO(2) holds index of column. C Internal variables CHARACTER INDFMT*16, MXTYPE*3, PTRFMT*16, VALFMI*20, VALFMO*20 LOGICAL YESA INTEGER DEC,I,IFLAG,INDN,INDCRD,J,K,NELTVL,NN,OUTFIL,PTRCRD, * PTRN,TOTCRD,VALCRD,ICT59(10),INFO59(10) EXTERNAL MC54B,MC54C,MC54D,MC54E,MC59A C Check for trivial error returns INFO(1) = 0 INFO(2) = 0 IF (M.LE.0) THEN INFO(1) = -1 INFO(2) = M GO TO 500 ENDIF IF (N.LE.0) THEN INFO(1) = -2 INFO(2) = N GO TO 500 ENDIF IF (NE.LE.0) THEN INFO(1) = -3 INFO(2) = NE GO TO 500 ENDIF IF (ICNTL(4).LT.0 .OR. ICNTL(4).GT.3) THEN INFO(1) = -4 INFO(2) = ICNTL(4) GO TO 500 ENDIF IF (ICNTL(5).LT.0 .OR. ICNTL(5).GT.2) THEN INFO(1) = -5 INFO(2) = ICNTL(5) GO TO 500 ENDIF IF (ICNTL(3).EQ.0 .AND. ICNTL(4).EQ.1) THEN INFO(1) = -6 GO TO 500 ENDIF IF (ICNTL(5).NE.1) THEN NN = N IF ((ICNTL(5).NE.1 .AND. ICNTL(5).NE.0) .AND. ICNTL(4).EQ.2) * NN = 2*N DO 10 I = 1,NN IF (IP(I+1).LT.IP(I)) THEN INFO(1) = -7 INFO(2) = I GO TO 500 ENDIF 10 CONTINUE ENDIF DO 20 I = 1,NE IF (IND(I).LE.0 .OR. IND(I).GT.M) THEN INFO(1) = -8 INFO(2) = I GO TO 500 ENDIF IF (ICNTL(5).EQ.1) THEN IF (IND(NE+I).LE.0 .OR. IND(NE+I).GT.N) THEN INFO(1) = -8 INFO(2) = I GO TO 500 ENDIF ENDIF 20 CONTINUE C If input is by columns, check that rows are in order. IF (ICNTL(5).EQ.0) THEN DO 40 J=1,N DO 30 K=IP(J),IP(J+1)-2 IF (IND(K).GE.IND(K+1)) THEN INFO(1) = -9 INFO(2) = J GO TO 500 ENDIF 30 CONTINUE 40 CONTINUE ENDIF C Set MXTYPE and determine form of input IF (ICNTL(3).EQ.0) THEN MXTYPE(1:1) = 'p' YESA = .FALSE. ELSE MXTYPE(1:1) = 'r' YESA = .TRUE. ENDIF IF (ICNTL(4).EQ.0) MXTYPE(2:2) = 's' IF (ICNTL(4).EQ.1) MXTYPE(2:2) = 'z' IF (ICNTL(4).EQ.2) MXTYPE(2:2) = 'r' IF (ICNTL(4).EQ.3) MXTYPE(2:2) = 'u' IF (ICNTL(5).EQ.0 .OR. ICNTL(5).EQ.1) THEN MXTYPE(3:3) = 'a' C Generate Col ptr/row index form if not input IF (ICNTL(5).EQ.1) THEN ICT59(1) = 1 ICT59(2) = 1 ICT59(3) = 1 IF (YESA) ICT59(3) = 0 ICT59(4) = -1 ICT59(5) = -1 ICT59(6) = 0 CALL MC59A(ICT59,N,M,NE,IND,NE,IND(NE+1),NE,VALUE, * MAX(M,N)+1,IP,MAX(M,N)+1,IW,INFO59) IFLAG = INFO59(1) ENDIF NELTVL = NE ELSE MXTYPE(3:3) = 'e' C Calculate NELTVL NELTVL = 0 DO 50 I=1,N IF (MXTYPE(2:2) .EQ. 'u') THEN NELTVL = NELTVL+(IP(I+1)-IP(I))*(IP(I+1)-IP(I)) ELSE IF (MXTYPE(2:2) .EQ. 'r') THEN C Rectangular element matrices NELTVL = NELTVL+((IP(2*I)-IP(2*I-1))*(IP(2*I+1)-IP(2*I))) ELSE C Symmetric or skew-symmetric element matrices NELTVL = NELTVL+((IP(I+1)-IP(I))*(IP(I+1)-IP(I)+1))/2 ENDIF 50 CONTINUE ENDIF OUTFIL = ICNTL(1) C Set format and number card images for pointer array CALL MC54D(NE+1,PTRN,PTRFMT) PTRCRD = N/PTRN + 1 IF (MXTYPE(3:3).EQ.'e' .AND. MXTYPE(2:2) .EQ. 'r') * PTRCRD = (2*N)/PTRN + 1 C Set format and number card images for index array CALL MC54D(M,INDN,INDFMT) INDCRD = (NE-1)/INDN + 1 C Set number of card images for matrix entries IF (YESA) THEN C DEC is set to number of significant decimal digits in matrix reals DEC = ICNTL(2) IF (DEC.LE.1 .OR. DEC.GT.17) DEC = 17 CALL MC54E(DEC,VALFMI,VALFMO,VALCRD) VALCRD = (NELTVL-1)/VALCRD + 1 ELSE VALCRD = 0 VALFMI = ' ' ENDIF TOTCRD = PTRCRD + INDCRD + VALCRD C Write header card IF (MXTYPE(3:3).EQ.'a') NELTVL = 0 WRITE( OUTFIL, 999 ) TITLE, KEY, 1 TOTCRD, PTRCRD, INDCRD, VALCRD, 2 MXTYPE, M, N, NE, NELTVL, 3 PTRFMT, INDFMT, VALFMI 999 FORMAT ( A72, A8 / I14, 3(1X,I13) / A3, 11X, 4(1X,I13) / 1 2A16, A20 ) C Write pointer array NN = N IF (MXTYPE(3:3).EQ.'e' .AND. MXTYPE(2:2) .EQ. 'r') NN = 2*N CALL MC54B(NN+1,IP,PTRFMT,OUTFIL) C Write indices CALL MC54B(NE,IND,INDFMT,OUTFIL) C Write matrix entries IF (MXTYPE(3:3).EQ.'a') NELTVL = NE IF (YESA) CALL MC54C(NELTVL,VALUE,VALFMO,OUTFIL) 500 RETURN END SUBROUTINE MC54B(N,IND,FMTIND,LUNIT) INTEGER N,LUNIT,IND(N) CHARACTER*16 FMTIND WRITE(LUNIT,FMTIND) IND RETURN END SUBROUTINE MC54C(N,A,FMTA,LUNIT) INTEGER N,LUNIT REAL A(N) CHARACTER*20 FMTA WRITE(LUNIT,FMTA) A RETURN END SUBROUTINE MC54D(N,NLIN,FMT) INTEGER N,NLIN,NN,K,NL(11) CHARACTER*16 FMT,FMTAB(11) DATA FMTAB / '(40I2)', '(26I3)', '(20I4)', '(16I5)', * '(13I6)', '(11I7)', '(10I8)', '(8I9)', * '(8I10)', '(7I11)', '(4I20)'/, * NL / 40,26,20,16,13,11,10,8,8,7,4 / C Note that any signed 32-bit integer can be held in a field of length 11, and C any signed 64-bit integer can be held in a field of length 20, both allowing C a space for the sign. NN = N DO 10 K=1,N IF (NN.LT.10) GO TO 20 NN = NN/10 10 CONTINUE 20 IF (K.LE.10) THEN FMT = FMTAB(K) NLIN = NL(K) ELSE FMT = FMTAB(11) NLIN = NL(11) ENDIF RETURN END SUBROUTINE MC54E(DEC,VALFMI,VALFMO,VALCRD) INTEGER DEC CHARACTER*20 VALFMI,VALFMO INTEGER VALCRD CHARACTER*20 FMT(16),FMT1(16) INTEGER LEN(16) DATA FMT / '(8E10.1E3)', '(7E11.2E3)', '(6E12.3E3)', * '(6E13.4E3)', * '(5E14.5E3)', '(5E15.6E3)', '(5E16.7E3)', * '(4E17.8E3)', '(4E18.9E3)', '(4E19.10E3)', * '(4E20.11E3)', '(3E21.12E3)', '(3E22.13E3)', * '(3E23.14E3)', '(3E24.15E3)', '(3E25.16E3)'/, * FMT1 / '(1P,8E10.1E3)', '(1P,7E11.2E3)', '(1P,6E12.3E3)', * '(1P,6E13.4E3)', * '(1P,5E14.5E3)', '(1P,5E15.6E3)', '(1P,5E16.7E3)', * '(1P,4E17.8E3)', '(1P,4E18.9E3)', '(1P,4E19.10E3)', * '(1P,4E20.11E3)', '(1P,3E21.12E3)', '(1P,3E22.13E3)', * '(1P,3E23.14E3)', '(1P,3E24.15E3)', '(1P,3E25.16E3)'/, * LEN / 8,7,6,6,5,5,5,4,4,4,4,3,3,3,3,3/ VALFMI = FMT(DEC-1) VALFMO = FMT1(DEC-1) VALCRD = LEN(DEC-1) RETURN END * COPYRIGHT (c) 2000 Council for the Central Laboratory * of the Research Councils *######DATE 22 Dec. 2000 SUBROUTINE MC59A(ICNTL,NC,NR,NE,IRN,LJCN,JCN,LA,A,LIP,IP, & LIW,IW,INFO) INTEGER LA,LIP,LIW,LJCN,NC,NE,NR REAL A(LA) INTEGER ICNTL(10),IP(LIP),INFO(10),IRN(NE),IW(LIW),JCN(LJCN) INTEGER I,ICNTL1,ICNTL2,ICNTL3,ICNTL6,LAA INTEGER IDUP,IOUT,IUP,JOUT,LP,MP,KNE,PART LOGICAL LCHECK EXTERNAL MC59B,MC59C,MC59D,MC59E,MC59F INTRINSIC MAX DO 10 I = 1,10 INFO(I) = 0 10 CONTINUE ICNTL1 = ICNTL(1) ICNTL2 = ICNTL(2) ICNTL3 = ICNTL(3) ICNTL6 = ICNTL(6) LCHECK = (ICNTL1.EQ.0) LP = ICNTL(4) MP = ICNTL(5) IF (ICNTL2.GT.2 .OR. ICNTL2.LT.0) THEN INFO(1) = -1 INFO(2) = ICNTL2 IF (LP.GT.0) THEN WRITE (LP,FMT=9000) INFO(1) WRITE (LP,FMT=9010) ICNTL2 END IF GO TO 70 END IF IF (ICNTL6.GT.2 .OR. ICNTL6.LT.-2) THEN INFO(1) = -11 INFO(2) = ICNTL6 IF (LP.GT.0) THEN WRITE (LP,FMT=9000) INFO(1) WRITE (LP,FMT=9150) ICNTL6 END IF GO TO 70 END IF IF (NC.LT.1) THEN INFO(1) = -2 INFO(2) = NC IF (LP.GT.0) THEN WRITE (LP,FMT=9000) INFO(1) WRITE (LP,FMT=9020) NC END IF GO TO 70 END IF IF (NR.LT.1) THEN INFO(1) = -3 INFO(2) = NR IF (LP.GT.0) THEN WRITE (LP,FMT=9000) INFO(1) WRITE (LP,FMT=9030) NR END IF GO TO 70 END IF IF (ICNTL6.NE.0 .AND. NR.NE.NC) THEN INFO(1) = -3 INFO(2) = NR IF (LP.GT.0) THEN WRITE (LP,FMT=9000) INFO(1) WRITE (LP,FMT=9035) NC,NR END IF GO TO 70 END IF IF (NE.LT.1) THEN INFO(1) = -4 INFO(2) = NE IF (LP.GT.0) THEN WRITE (LP,FMT=9000) INFO(1) WRITE (LP,FMT=9040) NE END IF GO TO 70 END IF IF (ICNTL2.EQ.0 .OR. ICNTL2.EQ.1) THEN IF (LJCN.LT.NE) THEN INFO(1) = -5 INFO(2) = NE END IF ELSE IF (LJCN.LT.1) THEN INFO(1) = -5 INFO(2) = 1 END IF END IF IF (INFO(1).EQ.-5) THEN IF (LP.GT.0) THEN WRITE (LP,FMT=9000) INFO(1) WRITE (LP,FMT=9050) LJCN,INFO(2) END IF GO TO 70 END IF IF (ICNTL3.EQ.0) THEN IF (LA.LT.NE) THEN INFO(1) = -6 INFO(2) = NE END IF ELSE IF (LA.LT.1) THEN INFO(1) = -6 INFO(2) = 1 END IF END IF IF (INFO(1).EQ.-6) THEN IF (LP.GT.0) THEN WRITE (LP,FMT=9000) INFO(1) WRITE (LP,FMT=9060) LA,INFO(2) END IF GO TO 70 END IF IF (ICNTL2.EQ.0 .OR. ICNTL2.EQ.2) THEN IF (LIP.LT.NC+1) THEN INFO(1) = -7 INFO(2) = NC+1 END IF ELSE IF (LIP.LT.MAX(NR,NC)+1) THEN INFO(1) = -7 INFO(2) = MAX(NR,NC)+1 END IF IF (INFO(1).EQ.-7) THEN IF (LP.GT.0) THEN WRITE (LP,FMT=9000) INFO(1) WRITE (LP,FMT=9065) LIP,INFO(2) END IF GO TO 70 END IF IF (LIW.LT.MAX(NR,NC)+1) THEN INFO(1) = -8 INFO(2) = MAX(NR,NC)+1 IF (LP.GT.0) THEN WRITE (LP,FMT=9000) INFO(1) WRITE (LP,FMT=9070) LIW,INFO(2) END IF GO TO 70 END IF LAA = NE IF (ICNTL3.NE.0) LAA = 1 IOUT = 0 JOUT = 0 IDUP = 0 IUP = 0 PART = 0 IF (ICNTL6.NE.0) PART = 1 IF (ICNTL2.EQ.0) THEN CALL MC59B(LCHECK,PART,NC,NR,NE,IRN,JCN,LAA,A,IP,IW, + IOUT,JOUT,KNE) IF (KNE.EQ.0) GO TO 50 IF (LCHECK) CALL MC59E(NC,NR,NE,IRN,LIP,IP,LAA,A,IW,IDUP, & KNE,ICNTL6) ELSE IF (ICNTL2.EQ.1) THEN IF (ICNTL6.NE.0) PART = -1 CALL MC59B(LCHECK,PART,NR,NC,NE,JCN,IRN,LAA,A,IW,IP, + JOUT,IOUT,KNE) IF (KNE.EQ.0) GO TO 50 IF (LCHECK) CALL MC59E(NR,NC,NE,JCN,NR+1,IW,LAA,A,IP, & IDUP,KNE,ICNTL6) CALL MC59C(NC,NR,KNE,IRN,JCN,LAA,A,IP,IW) ELSE IF (ICNTL2.EQ.2) THEN IF (LCHECK) THEN CALL MC59F(NC,NR,NE,IRN,NC+1,IP,LAA,A,LIW,IW,IDUP, + IOUT,IUP,KNE,ICNTL6,INFO) IF (INFO(1).EQ.-9) GO TO 40 IF (KNE.EQ.0) GO TO 50 ELSE KNE = NE END IF CALL MC59D(NC,KNE,IRN,IP,LAA,A) END IF INFO(3) = IDUP INFO(4) = IOUT INFO(5) = JOUT INFO(6) = KNE INFO(7) = IUP IF (IDUP.GT.0) INFO(1) = INFO(1) + 1 IF (IOUT.GT.0) INFO(1) = INFO(1) + 2 IF (JOUT.GT.0) INFO(1) = INFO(1) + 4 IF (INFO(1).GT.0 .AND. MP.GT.0) THEN WRITE (MP,FMT=9080) INFO(1) IF (IOUT.GT.0) WRITE (MP,FMT=9090) IOUT IF (JOUT.GT.0) WRITE (MP,FMT=9110) JOUT IF (IDUP.GT.0) WRITE (MP,FMT=9100) IDUP IF (IUP.GT.0) WRITE (MP,FMT=9130) IUP END IF GO TO 70 40 INFO(3) = IDUP INFO(4) = IOUT INFO(7) = IUP IF (LP.GT.0) THEN WRITE (LP,FMT=9000) INFO(1) WRITE (LP,FMT=9140) END IF GO TO 70 50 INFO(1) = -10 INFO(4) = IOUT INFO(5) = JOUT INFO(2) = IOUT + JOUT IF (LP.GT.0) THEN WRITE (LP,FMT=9000) INFO(1) WRITE (LP,FMT=9120) END IF 70 RETURN 9000 FORMAT (/,' *** Error return from MC59A *** INFO(1) = ',I3) 9010 FORMAT (1X,'ICNTL(2) = ',I2,' is out of range') 9020 FORMAT (1X,'NC = ',I6,' is out of range') 9030 FORMAT (1X,'NR = ',I6,' is out of range') 9035 FORMAT (1X,'Symmetric case. NC = ',I6,' but NR = ',I6) 9040 FORMAT (1X,'NE = ',I10,' is out of range') 9050 FORMAT (1X,'Increase LJCN from ',I10,' to at least ',I10) 9060 FORMAT (1X,'Increase LA from ',I10,' to at least ',I10) 9065 FORMAT (1X,'Increase LIP from ',I8,' to at least ',I10) 9070 FORMAT (1X,'Increase LIW from ',I8,' to at least ',I10) 9080 FORMAT (/,' *** Warning message from MC59A *** INFO(1) = ',I3) 9090 FORMAT (1X,I8,' entries in IRN supplied by the user were ', + /,' out of range and were ignored by the routine') 9100 FORMAT (1X,I8,' duplicate entries were supplied by the user') 9110 FORMAT (1X,I8,' entries in JCN supplied by the user were ', + /,' out of range and were ignored by the routine') 9120 FORMAT (1X,'All entries out of range') 9130 FORMAT (1X,I8,' of these entries were in the upper triangular ', + /,' part of matrix') 9140 FORMAT (1X,'Entries in IP are not monotonic increasing') 9150 FORMAT (1X,'ICNTL(6) = ',I2,' is out of range') END C*********************************************************************** SUBROUTINE MC59B(LCHECK,PART,NC,NR,NE,IRN,JCN,LA,A,IP,IW,IOUT, + JOUT,KNE) INTEGER LA,NC,NE,NR,IOUT,JOUT,KNE,PART LOGICAL LCHECK REAL A(LA) INTEGER IP(NC+1),IRN(NE),IW(NC+1),JCN(NE) REAL ACE,ACEP INTEGER I,ICE,ICEP,J,JCE,JCEP,K,L,LOC DO 10 J = 1,NC + 1 IW(J) = 0 10 CONTINUE KNE = 0 IOUT = 0 JOUT = 0 IF (LCHECK) THEN IF (LA.GT.1) THEN IF (PART.EQ.0) THEN DO 20 K = 1,NE I = IRN(K) J = JCN(K) IF (I.GT.NR .OR. I.LT.1) THEN IOUT = IOUT + 1 IF (J.GT.NC .OR. J.LT.1) JOUT = JOUT + 1 ELSE IF (J.GT.NC .OR. J.LT.1) THEN JOUT = JOUT + 1 ELSE KNE = KNE + 1 IRN(KNE) = I JCN(KNE) = J A(KNE) = A(K) IW(J) = IW(J) + 1 END IF 20 CONTINUE ELSE IF (PART.EQ.1) THEN DO 21 K = 1,NE I = IRN(K) J = JCN(K) IF (I.GT.NR .OR. I.LT.1) THEN IOUT = IOUT + 1 IF (J.GT.NC .OR. J.LT.1) JOUT = JOUT + 1 ELSE IF (J.GT.NC .OR. J.LT.1) THEN JOUT = JOUT + 1 ELSE KNE = KNE + 1 IF (I.LT.J) THEN IRN(KNE) = J JCN(KNE) = I IW(I) = IW(I) + 1 ELSE IRN(KNE) = I JCN(KNE) = J IW(J) = IW(J) + 1 END IF A(KNE) = A(K) END IF 21 CONTINUE ELSE IF (PART.EQ.-1) THEN DO 22 K = 1,NE I = IRN(K) J = JCN(K) IF (I.GT.NR .OR. I.LT.1) THEN IOUT = IOUT + 1 IF (J.GT.NC .OR. J.LT.1) JOUT = JOUT + 1 ELSE IF (J.GT.NC .OR. J.LT.1) THEN JOUT = JOUT + 1 ELSE KNE = KNE + 1 IF (I.GT.J) THEN IRN(KNE) = J JCN(KNE) = I IW(I) = IW(I) + 1 ELSE IRN(KNE) = I JCN(KNE) = J IW(J) = IW(J) + 1 END IF A(KNE) = A(K) END IF 22 CONTINUE END IF ELSE IF (PART.EQ.0) THEN DO 25 K = 1,NE I = IRN(K) J = JCN(K) IF (I.GT.NR .OR. I.LT.1) THEN IOUT = IOUT + 1 IF (J.GT.NC .OR. J.LT.1) JOUT = JOUT + 1 ELSE IF (J.GT.NC .OR. J.LT.1) THEN JOUT = JOUT + 1 ELSE KNE = KNE + 1 IRN(KNE) = I JCN(KNE) = J IW(J) = IW(J) + 1 END IF 25 CONTINUE ELSE IF (PART.EQ.1) THEN DO 26 K = 1,NE I = IRN(K) J = JCN(K) IF (I.GT.NR .OR. I.LT.1) THEN IOUT = IOUT + 1 IF (J.GT.NC .OR. J.LT.1) JOUT = JOUT + 1 ELSE IF (J.GT.NC .OR. J.LT.1) THEN JOUT = JOUT + 1 ELSE KNE = KNE + 1 IF (I.LT.J) THEN IRN(KNE) = J JCN(KNE) = I IW(I) = IW(I) + 1 ELSE IRN(KNE) = I JCN(KNE) = J IW(J) = IW(J) + 1 END IF END IF 26 CONTINUE ELSE IF (PART.EQ.-1) THEN DO 27 K = 1,NE I = IRN(K) J = JCN(K) IF (I.GT.NR .OR. I.LT.1) THEN IOUT = IOUT + 1 IF (J.GT.NC .OR. J.LT.1) JOUT = JOUT + 1 ELSE IF (J.GT.NC .OR. J.LT.1) THEN JOUT = JOUT + 1 ELSE KNE = KNE + 1 IF (I.GT.J) THEN IRN(KNE) = J JCN(KNE) = I IW(I) = IW(I) + 1 ELSE IRN(KNE) = I JCN(KNE) = J IW(J) = IW(J) + 1 END IF END IF 27 CONTINUE END IF END IF IF (KNE.EQ.0) GO TO 130 ELSE KNE = NE IF (PART.EQ.0) THEN DO 30 K = 1,NE J = JCN(K) IW(J) = IW(J) + 1 30 CONTINUE ELSE IF (PART.EQ.1) THEN DO 35 K = 1,NE I = IRN(K) J = JCN(K) IF (I.LT.J) THEN IRN(K) = J JCN(K) = I IW(I) = IW(I) + 1 ELSE IW(J) = IW(J) + 1 END IF 35 CONTINUE ELSE IF (PART.EQ.-1) THEN DO 36 K = 1,NE I = IRN(K) J = JCN(K) IF (I.GT.J) THEN IRN(K) = J JCN(K) = I IW(I) = IW(I) + 1 ELSE IW(J) = IW(J) + 1 END IF 36 CONTINUE END IF END IF IP(1) = 1 DO 37 J = 2,NC + 1 IP(J) = IW(J-1) + IP(J-1) IW(J-1) = IP(J-1) 37 CONTINUE IF (LA.EQ.1) THEN DO 70 L = 1,NC DO 60 K = IW(L),IP(L+1) - 1 ICE = IRN(K) JCE = JCN(K) DO 40 J = 1,NE IF (JCE.EQ.L) GO TO 50 LOC = IW(JCE) JCEP = JCN(LOC) ICEP = IRN(LOC) IW(JCE) = LOC + 1 JCN(LOC) = JCE IRN(LOC) = ICE JCE = JCEP ICE = ICEP 40 CONTINUE 50 JCN(K) = JCE IRN(K) = ICE 60 CONTINUE 70 CONTINUE ELSE DO 120 L = 1,NC DO 110 K = IW(L),IP(L+1) - 1 ICE = IRN(K) JCE = JCN(K) ACE = A(K) DO 90 J = 1,NE IF (JCE.EQ.L) GO TO 100 LOC = IW(JCE) JCEP = JCN(LOC) ICEP = IRN(LOC) IW(JCE) = LOC + 1 JCN(LOC) = JCE IRN(LOC) = ICE JCE = JCEP ICE = ICEP ACEP = A(LOC) A(LOC) = ACE ACE = ACEP 90 CONTINUE 100 JCN(K) = JCE IRN(K) = ICE A(K) = ACE 110 CONTINUE 120 CONTINUE END IF 130 CONTINUE RETURN END C********************************************************** SUBROUTINE MC59C(NC,NR,NE,IRN,JCN,LA,A,IP,IW) INTEGER LA,NC,NE,NR REAL A(LA) INTEGER IP(NC+1),IRN(NE),IW(NR+1),JCN(NE) REAL ACE,ACEP INTEGER I,ICE,ICEP,J,J1,J2,K,L,LOC,LOCP DO 10 J = 1,NC IP(J) = 0 10 CONTINUE IF (LA.GT.1) THEN DO 20 K = 1,NE I = JCN(K) IP(I) = IP(I) + 1 IRN(K) = JCN(K) 20 CONTINUE IP(NC+1) = NE + 1 IP(1) = IP(1) + 1 DO 30 J = 2,NC IP(J) = IP(J) + IP(J-1) 30 CONTINUE DO 50 I = NR,1,-1 J1 = IW(I) J2 = IW(I+1) - 1 DO 40 J = J1,J2 K = IRN(J) L = IP(K) - 1 JCN(J) = L IRN(J) = I IP(K) = L 40 CONTINUE 50 CONTINUE IP(NC+1) = NE + 1 DO 70 J = 1,NE LOC = JCN(J) IF (LOC.EQ.0) GO TO 70 ICE = IRN(J) ACE = A(J) JCN(J) = 0 DO 60 K = 1,NE LOCP = JCN(LOC) ICEP = IRN(LOC) ACEP = A(LOC) JCN(LOC) = 0 IRN(LOC) = ICE A(LOC) = ACE IF (LOCP.EQ.0) GO TO 70 ICE = ICEP ACE = ACEP LOC = LOCP 60 CONTINUE 70 CONTINUE ELSE DO 90 K = 1,NE I = JCN(K) IP(I) = IP(I) + 1 90 CONTINUE IP(NC+1) = NE + 1 IP(1) = IP(1) + 1 DO 100 J = 2,NC IP(J) = IP(J) + IP(J-1) 100 CONTINUE DO 120 I = NR,1,-1 J1 = IW(I) J2 = IW(I+1) - 1 DO 110 J = J1,J2 K = JCN(J) L = IP(K) - 1 IRN(L) = I IP(K) = L 110 CONTINUE 120 CONTINUE END IF RETURN END C********************************************************** SUBROUTINE MC59D(NC,NE,IRN,IP,LA,A) INTEGER LA,NC,NE REAL A(LA) INTEGER IRN(NE),IP(NC) REAL ACE INTEGER ICE,IK,J,JJ,K,KDUMMY,KLO,KMAX,KOR INTRINSIC ABS IF (LA.GT.1) THEN KMAX = NE DO 50 JJ = 1,NC J = NC + 1 - JJ KLO = IP(J) + 1 IF (KLO.GT.KMAX) GO TO 40 KOR = KMAX DO 30 KDUMMY = KLO,KMAX ACE = A(KOR-1) ICE = IRN(KOR-1) DO 10 K = KOR,KMAX IK = IRN(K) IF (ABS(ICE).LE.ABS(IK)) GO TO 20 IRN(K-1) = IK A(K-1) = A(K) 10 CONTINUE K = KMAX + 1 20 IRN(K-1) = ICE A(K-1) = ACE KOR = KOR - 1 30 CONTINUE 40 KMAX = KLO - 2 50 CONTINUE ELSE KMAX = NE DO 150 JJ = 1,NC J = NC + 1 - JJ KLO = IP(J) + 1 IF (KLO.GT.KMAX) GO TO 140 KOR = KMAX DO 130 KDUMMY = KLO,KMAX ICE = IRN(KOR-1) DO 110 K = KOR,KMAX IK = IRN(K) IF (ABS(ICE).LE.ABS(IK)) GO TO 120 IRN(K-1) = IK 110 CONTINUE K = KMAX + 1 120 IRN(K-1) = ICE KOR = KOR - 1 130 CONTINUE 140 KMAX = KLO - 2 150 CONTINUE END IF END C*********************************************************************** SUBROUTINE MC59E(NC,NR,NE,IRN,LIP,IP,LA,A,IW,IDUP,KNE,ICNTL6) INTEGER ICNTL6,IDUP,KNE,LIP,LA,NC,NR,NE REAL A(LA) INTEGER IRN(NE),IP(LIP),IW(NR) INTEGER I,J,K,KSTART,KSTOP,NZJ IDUP = 0 KNE = 0 DO 10 I = 1,NR IW(I) = 0 10 CONTINUE KSTART = IP(1) IF (LA.GT.1) THEN NZJ = 0 DO 30 J = 1,NC KSTOP = IP(J+1) IP(J+1) = IP(J) DO 20 K = KSTART,KSTOP - 1 I = IRN(K) IF (IW(I).LE.NZJ) THEN KNE = KNE + 1 IRN(KNE) = I A(KNE) = A(K) IP(J+1) = IP(J+1) + 1 IW(I) = KNE ELSE IDUP = IDUP + 1 IF (ICNTL6.GE.0) A(IW(I)) = A(IW(I)) + A(K) END IF 20 CONTINUE KSTART = KSTOP NZJ = KNE 30 CONTINUE ELSE DO 50 J = 1,NC KSTOP = IP(J+1) IP(J+1) = IP(J) DO 40 K = KSTART,KSTOP - 1 I = IRN(K) IF (IW(I).LT.J) THEN KNE = KNE + 1 IRN(KNE) = I IP(J+1) = IP(J+1) + 1 IW(I) = J ELSE IDUP = IDUP + 1 END IF 40 CONTINUE KSTART = KSTOP 50 CONTINUE END IF RETURN END C*********************************************************************** SUBROUTINE MC59F(NC,NR,NE,IRN,LIP,IP,LA,A,LIW,IW,IDUP,IOUT, + IUP,KNE,ICNTL6,INFO) INTEGER ICNTL6,IDUP,IOUT,IUP,KNE,LA,LIP,LIW,NC,NR,NE REAL A(LA) INTEGER IRN(NE),IP(LIP),IW(LIW),INFO(2) INTEGER I,J,K,KSTART,KSTOP,NZJ,LOWER IDUP = 0 IOUT = 0 IUP = 0 KNE = 0 DO 10 I = 1,NR IW(I) = 0 10 CONTINUE KSTART = IP(1) LOWER = 1 IF (LA.GT.1) THEN NZJ = 0 DO 30 J = 1,NC IF (ICNTL6.NE.0) LOWER = J KSTOP = IP(J+1) IF (KSTART.GT.KSTOP) THEN INFO(1) = -9 INFO(2) = J RETURN END IF IP(J+1) = IP(J) DO 20 K = KSTART,KSTOP - 1 I = IRN(K) IF (I.GT.NR .OR. I.LT.LOWER) THEN IOUT = IOUT + 1 IF (ICNTL6.NE.0 .AND. I.LT.J) IUP = IUP + 1 ELSE IF (IW(I).LE.NZJ) THEN KNE = KNE + 1 IRN(KNE) = I A(KNE) = A(K) IP(J+1) = IP(J+1) + 1 IW(I) = KNE ELSE IDUP = IDUP + 1 IF (ICNTL6.GE.0) A(IW(I)) = A(IW(I)) + A(K) END IF 20 CONTINUE KSTART = KSTOP NZJ = KNE 30 CONTINUE ELSE DO 50 J = 1,NC IF (ICNTL6.NE.0) LOWER = J KSTOP = IP(J+1) IF (KSTART.GT.KSTOP) THEN INFO(1) = -9 INFO(2) = J RETURN END IF IP(J+1) = IP(J) DO 40 K = KSTART,KSTOP - 1 I = IRN(K) IF (I.GT.NR .OR. I.LT.LOWER) THEN IOUT = IOUT + 1 IF (ICNTL6.NE.0 .AND. I.GT.1) IUP = IUP + 1 ELSE IF (IW(I).LT.J) THEN KNE = KNE + 1 IRN(KNE) = I IP(J+1) = IP(J+1) + 1 IW(I) = J ELSE IDUP = IDUP + 1 END IF 40 CONTINUE KSTART = KSTOP 50 CONTINUE END IF RETURN END