C COPYRIGHT (c) 1997 Council for the Central Laboratory * of the Research Councils C######DATE 3 September 1997 SUBROUTINE MC55AZ(ICNTL,TITLE,KEY,CASEID,DATTYP,POSITN,ORGNIZ, * M,NVEC,NAUXD,IP,IND,VALUE,INFO) INTEGER ICNTL(10) CHARACTER TITLE*72,KEY*8,DATTYP*3,POSITN*1,ORGNIZ*1, * CASEID*8 INTEGER M,NVEC,IP(*),IND(*) COMPLEX*16 VALUE(*) INTEGER NAUXD INTEGER INFO(5) C ================================================================ C Code for writing supplementary files in Rutherford-Boeing format. C For files with real entries MC55A/AD should be used. C For files with integer entries MC55AI 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 TITLE = descriptive title for matrix C KEY = identifier for associated matrix C CASEID = identifier for supplementary data C DATTYP = type of supplementary file C POSITN = position qualifier C ORGNIZ = data organization qualifier C M = row dimension of vector data C NVEC = number of vectors C NAUXD = total number of entries C IP = pointer array C IND = index array C VALUE = numerical values 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 unit 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 DATTYP is equal to ord, ipt, or icv, nor on a call to C MC55AI. C ICNTL(3) 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 supplementary data. This argument C is not altered 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 associated with this supplementary C data. This argument is not altered by the subroutine. C C CASEID is a CHARACTER*8 variable that must be set by the user to an C identifying key for the data. Different supplementary files of the C same type for the same matrix must have different values for CASEID. C This argument is not altered by the subroutine. C C DATTYP is a CHARACTER*3 variable that must be set by the user to C indicate the type of supplementary data being written. Current C possible values are C ord ... orderings C rhs ... right-hand sides C sln ... solutions C est ... estimates C evl ... eigenvalues C svl ... singular-values C evc ... eigenvectors C svc ... singular-vectors C sbv ... Schur basis vectors C sbm ... Schur basis matrix C sbp ... Schur basis parameters C ipt ... partition C icv ... covering C lvl ... Laplacian-values C lvc ... Laplacian-vectors C geo ... geometry C avl ... auxiliary values C C POSITN is a CHARACTER*1 that indicates position attributes of the C problem defined by the supplementary data. It need only be set for C values of DATTYP equal to ord, rhs, sln, est, evc, svc, lvc, ipt, icv, C and geo. C The meaning of POSITN in each of these cases is as follows: C ord .. POSITN has the value s, l, or r according to whether the C ordering is applied symmetrically, on the left, or on the C right. C rhs, sln, or est .. POSITN has the value l or r according to whether C the linear system is of the form Ax=b or A^Tx=b respectively. C evc, svc, lvc .. POSITN has the value s, l, or r according to whether C the vectors represent both left and right eigensolutions C (lvc and symmetric evc case), or the left or right vectors C respectively. C ipt, icv .. POSITN has the value s, l, or r according to whether the C partition or covering should be applied symmetrically, or C on the left (rows) or right (columns) respectively. C geo .. POSITN has the value s, l, or r according to whether the C geometry data applies symmetrically, or represents the C left/row vertices or right/column vertices respectively. C In these cases, this argument is not altered by the subroutine; C otherwise it is set to a blank character. C C ORGNIZ is a CHARACTER*1 variable that need only be set if DATTYP C is equal to rhs. It indicates whether the right-hand side is C sparse (s), dense (d), or elemental (e). In this case, this C argument is not altered by the subroutine; otherwise it is set C to a blank character. C C M is an INTEGER variable that must be set by the user to the row C dimension of the vector data. The actual meaning depends on the C value of DATTYP. For ord, rhs, sln, est, evc, svc, lvc, sbv, ipt, C icv, avl it is the row dimension of the vectors involved; for evl, C svl, lvl, sbm, sbp, geo it is the number of values, vectors, or C points involved. This argument is not altered by the subroutine. C C NVEC is an INTEGER variable that indicates the number of vectors C involved. It must be set by the user unless DATTYP equals evl, svl, C lvl, or sbp (in which case it is set by the subroutine to have the C value 1). In the case of DATTYP equal to geo, it should be set C to the geometric dimensionality. Unless set to 1, as above, this C argument is not altered by the subroutine. C C NAUXD is an INTEGER variable that need only be set if DATTYP is C equal to rhs and ORGNIZ is equal to e. In this case it must be set C to the total number of numerical entries in all right-hand sides; C In this case, this argument is not altered by the subroutine; in all C other cases it is altered. C C IP is an INTEGER array that need only be set by the user if C DATTYP is equal to rhs (ORGNIZ equal to s), ipt, or icv. For C rhs (ORGNIZ equal to s), it must be used to hold pointers to the start C of each sparse right-hand side in IND and VALUE; for ipt and icv C it must be set to point to the beginning of each partition or C covering in IND. C C IND is an INTEGER array that must be set by the user if DATTYP has C the value ord, rhs (ORGNIZ equal to s), ipt, or icv. For C ord, it must be set to the orderings; for C rhs (ORGNIZ equal to s), it must be set to the row indices of the C entries in the sparse right-hand sides; and for ipt and icv to the C indices for C each partitioning or covering. C C VALUE is a COMPLEX*16 array. It need not be set by the user if C DATTYP is equal to ord, ipt, or icv. Otherwise it must be set to C the value of the numerical data. C C INFO is an INTEGER array of length 5 that need not be set by the user. C On exit, INFO(1) will be zero if no error is detected, and negative C otherwise. Possible nonzero values are: C -1 M <= 0 C -2 NVEC <= 0 if DATTYP not equal to evl, svl, lvl, or sbp. C -3 Value of DATTYP not supported by this code C Internal variables C AUXFM1, AUXFM2, AUXFM3 are three CHARACTER*20 variables that are used C to hold formats for numerical data. For ord, rhs (with ORGNIZ equal C to d or e) C sln, est, evl, svl, lvl, evc, svc, lvc, sbv, sbm, sbp, geo, avl only C AUXFM1 is used and it must be set to the format for the numerical C entries. For DATTYP equal to rhs (ORGNIZ equal to s), AUXFM1 is set C to the format for right-hand side column pointers, AUXFM2 to the C format C for row indices of right-hand-sides and AUXFM3 to the format for C the numerical values of the right-hand sides; for ipt or ipv, C AUXFM1 must be set to the format for integer values in the pointer C array C and AUXFM2 to the format for the integer values in the index list C array. C CHARACTER NUMERF*1 CHARACTER*20 AUXFM1,AUXFM2,AUXFM3,FM1,FM3 INTEGER DEC,OUTFIL EXTERNAL MC55BZ,MC55CZ,MC55DZ,MC55EZ INFO(1) = 0 INFO(2) = 0 C Check for trivial error returns IF (M.LE.0) THEN INFO(1) = -1 INFO(2) = M GO TO 500 ENDIF C Set value of NVEC IF (DATTYP.EQ.'evl' .OR. DATTYP.EQ.'svl' .OR. 1 DATTYP.EQ.'lvl' .OR. DATTYP.EQ.'sbp') NVEC =1 C Test value of NVEC IF (NVEC.LE.0) THEN INFO(1) = -2 INFO(2) = NVEC GO TO 500 ENDIF IF (DATTYP.NE.'ord' .AND. DATTYP.NE.'rhs' .AND. 1 DATTYP.NE.'sln' .AND. DATTYP.NE.'est' .AND. 2 DATTYP.NE.'evl' .AND. DATTYP.NE.'svl' .AND. 3 DATTYP.NE.'evc' .AND. DATTYP.NE.'svc' .AND. 4 DATTYP.NE.'sbv' .AND. DATTYP.NE.'sbm' .AND. 5 DATTYP.NE.'sbp' .AND. DATTYP.NE.'ipt' .AND. 6 DATTYP.NE.'icv' .AND. DATTYP.NE.'lvl' .AND. 7 DATTYP.NE.'lvc' .AND. DATTYP.NE.'geo' .AND. 8 DATTYP.NE.'avl') THEN INFO(1) = -3 GO TO 500 ENDIF OUTFIL = ICNTL(1) IF (DATTYP.EQ.'evl' .OR. DATTYP.EQ.'svl' .OR. 1 DATTYP.EQ.'lvl' .OR. DATTYP.EQ.'sbv' .OR. 2 DATTYP.EQ.'sbm' .OR. DATTYP.EQ.'sbp' .OR. 3 DATTYP.EQ.'avl') POSITN = ' ' IF (DATTYP.NE.'rhs') ORGNIZ = ' ' C Note that, for INTEGER or REAL numerical data, MC55AI or C MC55A/AD should be used, respectively NUMERF = 'c' IF (DATTYP.EQ.'ord') NUMERF = 'i' IF ((DATTYP.EQ.'ipt') .OR. (DATTYP.EQ.'icv')) NUMERF = 'p' IF (ORGNIZ.NE.'e') THEN NAUXD = 0 IF (ORGNIZ .EQ. 's' .OR. DATTYP .EQ. 'icv' .OR. 1 DATTYP .EQ. 'ipt') NAUXD = IP(NVEC+1)-1 IF (ORGNIZ .EQ.'d') NAUXD = M*NVEC ENDIF C Set format statements AUXFM1 = ' ' AUXFM2 = ' ' AUXFM3 = ' ' FM1 = ' ' FM3 = ' ' IF (DATTYP .EQ. 'ord' .OR. DATTYP .EQ. 'ipt' .OR. * DATTYP .EQ. 'icv') THEN IF (DATTYP .EQ. 'ipt' .OR. DATTYP .EQ. 'icv') THEN CALL MC55DZ(NAUXD+1,AUXFM1) CALL MC55DZ(M,AUXFM2) ELSE CALL MC55DZ(M,AUXFM1) ENDIF ELSE DEC = ICNTL(2) IF (DEC.LT.2 .OR. DEC.GT.17) DEC = 17 IF (DATTYP .EQ. 'rhs' .AND. ORGNIZ.EQ.'s') THEN CALL MC55DZ(NAUXD+1,AUXFM1) CALL MC55DZ(M,AUXFM2) CALL MC55EZ(DEC,AUXFM3,FM3) ELSE CALL MC55EZ(DEC,AUXFM1,FM1) ENDIF ENDIF C Write header cards WRITE( OUTFIL, 999 ) TITLE, KEY, 1 DATTYP, POSITN, ORGNIZ, CASEID, 2 NUMERF, M, NVEC, NAUXD, 3 AUXFM1, AUXFM2, AUXFM3 999 FORMAT ( A72, A8 / A3, 2A1, 1X, A8, 1X, A1, 3(1X, I13) / 1 3A20 ) IF ((DATTYP.EQ.'rhs' .AND. ORGNIZ.EQ.'s') .OR. 1 (DATTYP.EQ.'ipt') .OR. (DATTYP.EQ.'icv')) THEN C Write pointer array CALL MC55BZ(NVEC+1,IP,AUXFM1,OUTFIL) C Write indices CALL MC55BZ(NAUXD,IND,AUXFM2,OUTFIL) IF (DATTYP.EQ.'ipt' .OR. DATTYP.EQ.'icv') GO TO 500 ENDIF IF (DATTYP .NE. 'rhs') NAUXD = M*NVEC IF (DATTYP.EQ.'ord') THEN CALL MC55BZ(NAUXD,IND,AUXFM1,OUTFIL) GO TO 500 ENDIF C Write values IF (DATTYP .EQ. 'rhs' .AND. ORGNIZ .EQ. 's') THEN CALL MC55CZ(NAUXD,VALUE,FM3,OUTFIL) ELSE CALL MC55CZ(NAUXD,VALUE,FM1,OUTFIL) ENDIF 500 RETURN END SUBROUTINE MC55BZ(N,IND,FMTIND,LUNIT) INTEGER N,LUNIT,IND(N) CHARACTER*20 FMTIND WRITE(LUNIT,FMTIND) IND RETURN END SUBROUTINE MC55CZ(N,A,FMTA,LUNIT) INTEGER N,LUNIT COMPLEX*16 A(N) CHARACTER*20 FMTA WRITE(LUNIT,FMTA) A RETURN END SUBROUTINE MC55DZ(N,FMT) INTEGER N,NN,K CHARACTER*16 FMT,FMTAB(12) DATA FMTAB / '(40I2)', '(26I3)', '(20I4)', '(16I5)', * '(13I6)', '(11I7)', '(10I8)', '(8I9)', * '(8I10)', '(7I11)', '(6I12)', '(4I20)'/ 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.11) THEN FMT = FMTAB(K) ELSE FMT = FMTAB(12) ENDIF RETURN END SUBROUTINE MC55EZ(DEC,VALFMI,VALFMO) INTEGER DEC CHARACTER*20 VALFMI,VALFMO CHARACTER*20 FMT(16),FMT1(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)'/ VALFMI = FMT(DEC-1) VALFMO = FMT1(DEC-1) RETURN END