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 MC54AZ(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(*)
COMPLEX*16 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 real entries MC54A/AD 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), unsymmetric (ICNTL(4)=3),
C or Hermitian (ICNTL(4)=4).
C A value of ICNTL(4)=1 is invalid if ICNTL(3)=0.
C An error flag is set if ICNTL(4) < 0 or > 4.
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 COMPLEX*16 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 MC54BZ,MC54CZ,MC54DZ,MC54EZ,MC59AZ
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.4) 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) = 'c'
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(4).EQ.4) MXTYPE(2:2) = 'h'
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 MC59AZ(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, skew, or Hermitian 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 MC54DZ(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 MC54DZ(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 MC54EZ(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 MC54BZ(NN+1,IP,PTRFMT,OUTFIL)
C Write indices
CALL MC54BZ(NE,IND,INDFMT,OUTFIL)
C Write matrix entries
IF (MXTYPE(3:3).EQ.'a') NELTVL = NE
IF (YESA) CALL MC54CZ(NELTVL,VALUE,VALFMO,OUTFIL)
500 RETURN
END
SUBROUTINE MC54BZ(N,IND,FMTIND,LUNIT)
INTEGER N,LUNIT,IND(N)
CHARACTER*16 FMTIND
WRITE(LUNIT,FMTIND) IND
RETURN
END
SUBROUTINE MC54CZ(N,A,FMTA,LUNIT)
INTEGER N,LUNIT
COMPLEX*16 A(N)
CHARACTER*20 FMTA
WRITE(LUNIT,FMTA) A
RETURN
END
SUBROUTINE MC54DZ(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 MC54EZ(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)', '(6E11.2E3)', '(6E12.3E3)',
* '(6E13.4E3)',
* '(4E14.5E3)', '(4E15.6E3)', '(4E16.7E3)',
* '(4E17.8E3)', '(4E18.9E3)', '(4E19.10E3)',
* '(4E20.11E3)', '(2E21.12E3)', '(2E22.13E3)',
* '(2E23.14E3)', '(2E24.15E3)', '(2E25.16E3)'/,
* FMT1 / '(1P,8E10.1E3)', '(1P,6E11.2E3)', '(1P,6E12.3E3)',
* '(1P,6E13.4E3)',
* '(1P,4E14.5E3)', '(1P,4E15.6E3)', '(1P,4E16.7E3)',
* '(1P,4E17.8E3)', '(1P,4E18.9E3)', '(1P,4E19.10E3)',
* '(1P,4E20.11E3)', '(1P,2E21.12E3)', '(1P,2E22.13E3)',
* '(1P,2E23.14E3)', '(1P,2E24.15E3)', '(1P,2E25.16E3)'/,
* LEN / 4,3,3,3,2,2,2,2,2,2,2,1,1,1,1,1/
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 MC59AZ(ICNTL,NC,NR,NE,IRN,LJCN,JCN,LA,A,LIP,IP,
& LIW,IW,INFO)
INTEGER LA,LIP,LIW,LJCN,NC,NE,NR
COMPLEX*16 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 MC59BZ,MC59CZ,MC59DZ,MC59EZ,MC59FZ
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 (ICNTL2.EQ.0) THEN
IF (ABS(ICNTL6).EQ.1) PART = 1
IF (ABS(ICNTL6).EQ.2) PART = 2
CALL MC59BZ(LCHECK,PART,NC,NR,NE,IRN,JCN,LAA,A,IP,IW,
+ IOUT,JOUT,KNE)
IF (KNE.EQ.0) GO TO 50
IF (LCHECK) CALL MC59EZ(NC,NR,NE,IRN,LIP,IP,LAA,A,IW,IDUP,
& KNE,ICNTL6)
ELSE IF (ICNTL2.EQ.1) THEN
IF (ABS(ICNTL6).EQ.1) PART = -1
IF (ABS(ICNTL6).EQ.2) PART = -2
CALL MC59BZ(LCHECK,PART,NR,NC,NE,JCN,IRN,LAA,A,IW,IP,
+ JOUT,IOUT,KNE)
IF (KNE.EQ.0) GO TO 50
IF (LCHECK) CALL MC59EZ(NR,NC,NE,JCN,NR+1,IW,LAA,A,IP,
& IDUP,KNE,ICNTL6)
CALL MC59CZ(NC,NR,KNE,IRN,JCN,LAA,A,IP,IW)
ELSE IF (ICNTL2.EQ.2) THEN
IF (LCHECK) THEN
CALL MC59FZ(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 MC59DZ(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 MC59AZ *** 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 MC59AZ *** 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 MC59BZ(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
COMPLEX*16 A(LA)
INTEGER IP(NC+1),IRN(NE),IW(NC+1),JCN(NE)
COMPLEX*16 ACE,ACEP
INTEGER I,ICE,ICEP,J,JCE,JCEP,K,L,LOC
INTRINSIC CONJG
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.2) THEN
DO 211 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
A(KNE) = CONJG(A(K))
ELSE
IRN(KNE) = I
JCN(KNE) = J
IW(J) = IW(J) + 1
A(KNE) = A(K)
END IF
END IF
211 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
ELSE IF (PART.EQ.-2) THEN
DO 221 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
write (6,*) 'swap'
IRN(KNE) = J
JCN(KNE) = I
IW(I) = IW(I) + 1
A(KNE) = CONJG(A(K))
ELSE
IRN(KNE) = I
JCN(KNE) = J
IW(J) = IW(J) + 1
A(KNE) = A(K)
END IF
END IF
221 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 .OR. PART.EQ.2) 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 .OR. PART.EQ.-2) 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 .OR. PART.EQ.2) 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
ELSE IF (PART.EQ.-2) THEN
IF (LA.GT.1) THEN
DO 37 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
A(K) = CONJG(A(K))
ELSE
IW(J) = IW(J) + 1
END IF
37 CONTINUE
ELSE
DO 38 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
38 CONTINUE
END IF
END IF
END IF
IP(1) = 1
DO 39 J = 2,NC + 1
IP(J) = IW(J-1) + IP(J-1)
IW(J-1) = IP(J-1)
39 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 MC59CZ(NC,NR,NE,IRN,JCN,LA,A,IP,IW)
INTEGER LA,NC,NE,NR
COMPLEX*16 A(LA)
INTEGER IP(NC+1),IRN(NE),IW(NR+1),JCN(NE)
COMPLEX*16 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 MC59DZ(NC,NE,IRN,IP,LA,A)
INTEGER LA,NC,NE
COMPLEX*16 A(LA)
INTEGER IRN(NE),IP(NC)
COMPLEX*16 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 MC59EZ(NC,NR,NE,IRN,LIP,IP,LA,A,IW,IDUP,KNE,ICNTL6)
INTEGER ICNTL6,IDUP,KNE,LIP,LA,NC,NR,NE
COMPLEX*16 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 MC59FZ(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
COMPLEX*16 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