C INCLUDE 'CALPGM.INC' C miscellaneous subroutines: ULIB.FOR C C Copyright (C) 1989, California Institute of Technology C All rights reserved. U. S. Government Sponsorship under C NASA Contract NAS7-918 is acknowledged. C C Herbert M. Pickett, 20 March 1989 C SUBROUTINE ORDBLK(NDM,NN,T,E,NBLK,ISBLK,P,IPASGN,IORD) INTEGER NDM,NN,NBLK,IORD INTEGER*2 ISBLK(*),IPASGN(*) REAL*8 T(NDM,*),E(*),P(*) C C SUBROUTINE TO ARRANGE EIGENVECTORS AND EIGENVALUES SO THE C BLOCK PROJECTION MATRIX IS DOMINANT DIAGONAL C AND ENERGIES ARE ORDERED WITHIN BLOCK C ON INPUT - C C T IS THE EIGENVECTOR MATRIX AND E IS THE EIGENVALUE C NDM IS THE DIMENSIONED LENGTH OF T AND NN IS THE ACTUAL LENGTH C NBLK IS THE NUMBER OF BLOCKS C ISBLK CONTAINS THE STARTING INDEX OF A BLOCK C (NOTE: ISBLK(NBLK+1)=NN+1 ) C C IORD <0 MEADS REVERSE SORT , IORD = 0 MEANS NO SORT C C RETURNS - C C SORTED T AND E C C P IS THE PROJECTION OF THE STATE ON ITS BLOCK C C IPASGN IS A SCRATCH VECTOR WHICH GIVES THE BLOCK ASSIGNMENT C POSITIVE VALUE MEANS PROVISIONAL ASSIGNMENT C ZERO VALUE MEANS FINAL ASSIGNMENT C C INTEGER I1,EIG,N,IZ,IBK,IS,NS,NLEFT,EIGMAX,ISWAP,IBIG,K REAL*8 DDOT,ETMP,PTMP,PCMP,PBIG,HALF PARAMETER (I1=1,HALF=0.5D0) C IF(NN.LE.1) THEN P(1)=1 RETURN ENDIF N=NN C C FIND LARGEST PROJECTION FOR EACH EIGENVECTOR C IZ=NBLK-1 C SKIP PROJECTION FINDING FOR ONLY 1 BLOCK AND GO TO ENERGY ORDERING IF(IZ.GT.0) THEN DO 3 EIG=1,N PBIG=1 PCMP=0 DO 4 IBK=1,IZ IS=ISBLK(IBK) NS=ISBLK(IBK+1)-IS PTMP=DDOT(NS,T(IS,EIG),I1,T(IS,EIG),I1) C .. ACCUMULATE PROJECTION FOR LAST BLOCK SPECIALLY PBIG=PBIG-PTMP IF(PTMP.GT.PCMP) THEN IPASGN(EIG)=IBK PCMP=PTMP C ..CHECK FOR AUTOMATIC LARGEST PROJECTION IF(PCMP.GE.HALF) GO TO 3 ENDIF 4 CONTINUE IF(PBIG.GT.PCMP) THEN IPASGN(EIG)=NBLK PCMP=PBIG ENDIF 3 P(EIG)=PCMP C C SEARCH FOR UNASSIGNED STATES C NLEFT=N-1 C DO WHILE(NLEFT.GT.0) 10 IF(NLEFT.GT.0) THEN PCMP=0 C ...FIND NEXT UNASSIGNED STATE WITH LARGEST PROJECTION DO 11 EIG=1,N IF(IPASGN(EIG).NE.0) THEN IF(P(EIG).GT.PCMP) THEN PCMP=P(EIG) EIGMAX=EIG ENDIF ENDIF 11 CONTINUE IBK=IPASGN(EIGMAX) IS=ISBLK(IBK) IF(IPASGN(IS).EQ.0) THEN C .. NO PLACE FOUND IN BLOCK SO FIND NEXT BEST PROJECTION PCMP=0 DO 14 IBK=1,NBLK IS=ISBLK(IBK) IF(IPASGN(IS).NE.0) THEN NS=ISBLK(IBK+1)-IS PTMP=DDOT(NS,T(IS,EIGMAX),I1,T(IS,EIGMAX),I1) IF(PTMP.GT.PCMP) THEN IPASGN(EIGMAX)=IBK PCMP=PTMP ENDIF ENDIF 14 CONTINUE P(EIGMAX)=PCMP ELSE IZ=ISBLK(IBK+1)-1 C .. FIND NEXT AVAILABLE POSITION IN BLOCK (PREFER PRESENT) C .. LAST AVAILABLE POSITION IS FIRST BLOCK POSITION IF(EIGMAX.GT.IS.AND.EIGMAX.LE.IZ) IZ=EIGMAX DO 12 ISWAP=IZ,IS,-1 IF(IPASGN(ISWAP).NE.0) GO TO 13 12 CONTINUE 13 IF(ISWAP.NE.EIGMAX) THEN IPASGN(EIGMAX)=IPASGN(ISWAP) PTMP=P(EIGMAX) P(EIGMAX)=P(ISWAP) P(ISWAP)=PTMP ETMP=E(EIGMAX) E(EIGMAX)=E(ISWAP) E(ISWAP)=ETMP CALL DSWAP(N,T(1,EIGMAX),I1,T(1,ISWAP),I1) ENDIF IPASGN(ISWAP)=0 NLEFT=NLEFT-1 ENDIF GO TO 10 ENDIF ENDIF C C ENERGY ORDER BLOCKS BY MAXIMUM SEARCH AND SWAP C IF(IORD.NE.0) THEN IF(IORD.LT.0) THEN DO 49 EIG=1,N 49 E(EIG)=-E(EIG) ENDIF DO 50 IBK=1,NBLK IS=ISBLK(IBK) IZ=ISBLK(IBK+1)-2 DO 55 EIG=IZ,IS,-1 IBIG=0 ISWAP=EIG+1 ETMP=E(ISWAP) C .. FIND MAX E DO 53 K=IS,EIG IF(E(K).GT.ETMP) THEN IBIG=K ETMP=E(K) ENDIF 53 CONTINUE IF(IBIG.NE.0) THEN C ...PERFORM SWAP WITH MAX E(IBIG)=E(ISWAP) E(ISWAP)=ETMP PTMP=P(IBIG) P(IBIG)=P(ISWAP) P(ISWAP)=PTMP CALL DSWAP(N,T(1,IBIG),I1,T(1,ISWAP),I1) ENDIF 55 CONTINUE 50 CONTINUE IF(IORD.LT.0) THEN DO 60 EIG=1,N 60 E(EIG)=-E(EIG) ENDIF ENDIF C set phase for positive diagonal PTMP=-1 DO 101 EIG=1,N IF( T(EIG,EIG).LT.0 ) CALL DSCAL(N,PTMP,T(1,EIG),I1) 101 CONTINUE RETURN END SUBROUTINE HDIAG(NM,NX,Z,D,E,IERR) C C INTEGER NM,NX,IERR REAL*8 Z(NM,*),D(*),E(*) C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED2, C NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). C C MODIFIED TO SPEED UP FOR SPARCE MATRICIES C C THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX TO A C SYMMETRIC TRIDIAGONAL MATRIX USING AND ACCUMULATING C ORTHOGONAL SIMILARITY TRANSFORMATIONS. C C THE SUBROUTINE THEN CALLS TRIAG WHICH DIAGONALIZES C C ON INPUT- C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT, C C NX IS THE ORDER OF THE MATRIX, C C Z CONTAINS THE REAL SYMMETRIC INPUT MATRIX. ONLY THE C LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED. C C ON CALL TO TRIAG - C C D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX, C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL C MATRIX IN ITS FIRST N-1 POSITIONS. E(N) IS ARBITRARY, C C Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX C PRODUCED IN THE REDUCTION. C C ON OUTPUT- C C D CONTAINS THE EIGENVALUES IN ARBITRARY ORDER C C Z CONTAINS THE EIGENVECTORS C C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE J-TH EIGENVALUE HAS NOT BEEN C DETERMINED AFTER 30 ITERATIONS. C C ------------------------------------------------------------------ C INTEGER I0,I1 PARAMETER (I0=0,I1=1) REAL*8 MACHEP,F,G,H,T,BB,DDOT,ONE,ZERO(2) PARAMETER (ONE=1.) INTEGER I,J,K,L,IZ,LZ,NT1,NT2,NDM,N SAVE ZERO DATA ZERO /0.,0./ MACHEP=1.E-30 N=NX NDM=NM C ********** FOR I=N STEP -1 UNTIL 3 DO -- ********** DO 300 I = N,3,-1 L = I - 1 J = L - 1 CALL DCOPY(L,Z(I,1),NDM,D,I1) F=D(L) IZ=0 H=0 C GET SUM OF SQUARES OF ELEMENTS AND FIND FIRST NON-ZERO ELEMENT DO 220 K=J,1,-1 G=D(K)*D(K) IF(G.GT.MACHEP) THEN IZ=K H=H+G ELSE D(K)=0 ENDIF 220 CONTINUE E(I)=F IF(IZ.NE.0) THEN T=F*F IF(H.LT.MACHEP*T) IZ=0 ENDIF IF(IZ.EQ.0) THEN D(L)=0 ELSE G= SIGN( SQRT( H +T ), F) E(I)= -G F= F + G D(L)=F BB=SIGN( ONE, F ) / SQRT( F*G ) LZ=I-IZ CALL DSCAL(LZ,BB,D(IZ),I1) CALL DCOPY(L,D,I1,Z(1,I),I1) C DO 230 J= 1, IZ 230 E(J)= DDOT(LZ,Z(IZ,J),I1,D(IZ),I1) F=-E(IZ)*D(IZ) DO 240 J = IZ+1, L NT1=J-IZ NT2=I-J T= DDOT(NT1,Z(J,IZ),NDM,D(IZ),I1) + +DDOT(NT2,Z(J, J), I1,D( J),I1) E(J)=T 240 F=F-T*D(J) C C ********** FORM REDUCED A ********** DO 260 J = IZ, L H = -D(J) G = F * H - E(J) K = I-J IF(ABS(H).GT.MACHEP) CALL DAXPY(K,H,E(J),I1,Z(J,J),I1) CALL DAXPY(K,G,D(J),I1,Z(J,J),I1) 260 CONTINUE DO 265 J = 1,IZ-1 G = -E(J) CALL DAXPY(LZ,G,D(IZ),I1,Z(IZ,J),I1) 265 CONTINUE ENDIF 300 CONTINUE C D(1)=Z(1,1) Z(1,1)=ONE IF(N.LE.1) THEN IERR=0 RETURN ENDIF D(N)=0 E(2)=Z(2,1) C ********** ACCUMULATION OF TRANSFORMATION MATRICES ********** DO 500 L = 2, N I = L + 1 J = L - 1 E(J)=E(L) F=D(L) D(L)=Z(L,L) IF(F.GT.MACHEP) THEN DO 360 K=1,J G= -DDOT(J,Z(1,I),I1,Z(1,K),I1) IF(ABS(G).GT.MACHEP) CALL DAXPY(J,G,Z(1,I),I1 + ,Z(1,K),I1) Z(L,K)= F * G Z(K,L)= - F * Z(K,I) 360 CONTINUE Z(L,L) = ONE -F*F ELSE Z(L,L) = ONE CALL DCOPY(J,ZERO,I0,Z(1,L),I1) CALL DCOPY(J,ZERO,I0,Z(L,1),NDM) ENDIF 500 CONTINUE C CALL TRIAG(NM,N,N,Z,D,E,IERR) RETURN END SUBROUTINE TRIAG(NM,NX,NZ,Z,D,E,IERR) C INTEGER NM,NX,NZ,IERR REAL*8 Z(0:*),D(*),E(*) C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE imtql2, C NUM. MATH. 12, 377-383(1968) BY MARTIN AND WILKINSON. C C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS C OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE QL METHOD. C C ON INPUT- C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT, C C NX IS THE ORDER OF THE MATRIX, C C NZ IS THE ACTUAL COLUMN LENGTH OF Z, C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX, C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX C IN ITS FIRST N-1 POSITIONS. E(N) IS ARBITRARY, C C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE C REDUCTION BY TRED2, IF PERFORMED. IF THE EIGENVECTORS C OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN C THE IDENTITY MATRIX. C C ON OUTPUT- C C D CONTAINS THE EIGENVALUES IN ARBITRARY ORDER. IF AN C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT C UNORDERED FOR INDICES 1,2,...,IERR-1, C C E HAS BEEN DESTROYED, C C Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC C TRIDIAGONAL (OR FULL) MATRIX. IF AN ERROR EXIT IS MADE, C Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED C EIGENVALUES, C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE J-TH EIGENVALUE HAS NOT BEEN C DETERMINED AFTER 30 ITERATIONS. C C ------------------------------------------------------------------ C C ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING C THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC. C C ********** INTEGER I1 PARAMETER (I1=1) REAL*8 MACHEP,F,B,P,G,R,C,S INTEGER*4 INDX1,INDX2 INTEGER N,NN,ITR,M,I,NDM,L,MM MACHEP=1.E-15 N=NX C IERR = 0 IF (N .LE. 1) RETURN NN=N-1 NDM=NM E(N) = 0 F = 0 B = 0 C DO 240 L = 1, NN ITR = 31 C ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT ********** DO WHILE(ITR.NE.0) M=N DO 110 I = L, NN IF (ABS(E(I)).LE.MACHEP*(ABS(D(I))+ABS(D(I+1)))) THEN M=I GO TO 120 ENDIF 110 CONTINUE C 120 P=D(L) IF (M .EQ. L) GO TO 240 ITR = ITR - 1 IF (ITR .EQ. 0) THEN C ** SET ERROR -- NO CONVERGENCE TO AN EIGENVALUE AFTER 30 ITERATIONS ** IERR = L RETURN ENDIF C ********** FORM SHIFT ********** G = E(L) G = ( D(L+1)-P )/( G + G ) R = SQRT(1+G*G) G = D(M) - P + E(L) / ( G +SIGN(R,G) ) C ********** QL TRANSFORMATION ********** P = D(M) C = 1 S = C C ********** FOR I=M-1 STEP -1 UNTIL L DO -- ********** MM=M-1 INDX2=NDM*MM DO 200 I = MM,L,-1 B=E(I)*C F=E(I)*S IF(ABS(F).GE.ABS(G)) THEN C=G/F R=SQRT(1+C*C) E(I+1)=F*R S=1/R C=C*S ELSE S=F/G R=SQRT(1+S*S) E(I+1)=G*R C=1/R S=S*C ENDIF F=C*D(I)-S*B G=C*B-S*P R=D(I)+P P=C*F-S*G G=S*F+C*G D(I+1)=R-P C ********** FORM VECTOR ********** INDX1=INDX2 INDX2=INDX2-NDM CALL DROT(NZ,Z(INDX1),I1,Z(INDX2),I1,C,S) C 200 CONTINUE C D(L) = P E(L) = G E(M) = 0 ENDDO 240 CONTINUE RETURN END C FUNCTION GETPAR(LU,NTOP,VAR,ERPAR,NPAR,IDPAR,PAR,LUBIN,LUOUT,PLBL) C C SUBROUTINE TO READ PARAMETERS AND VARIANCE C INTEGER GETPAR,LU,NPAR,LUBIN,LUOUT,IERR C LUBIN=0 for CALFIT INTEGER*4 NTOP,NL REAL*8 IDPAR(0:*),VAR(0:*),ERPAR(*),PAR(*), : ERMIN,VEC(3),ERMAX,VAL,ZERO(2) CHARACTER*10 PLBL(*),TLBL CHARACTER CARD*80 INTEGER I,N,I0,I1,IC,NPARTOT,GETVAR,NVAL,PCARD PARAMETER (I0=0,I1=1) SAVE ZERO DATA ZERO/0.,0./ GETPAR=0 ERMIN=1.E-19 ERMAX=ERMIN IF(LUBIN.EQ.0) ERMAX=1.E+38 N=NPAR NPAR=0 NPARTOT=0 DO 10 I=1,N VEC(1)=0 VEC(2)=0 VEC(3)=ERMAX READ (LU,100,IOSTAT=IERR) CARD IF(IERR.NE.0) GO TO 50 NVAL=3 IF(PCARD(CARD,VEC,NVAL).EQ.0) GO TO 50 IF(I.EQ.1) VEC(1)=ABS( VEC(1) ) IDPAR(I)=VEC(1) IF(VEC(1).GE.0) NPAR=NPAR+1 PAR(I)=VEC(2) ERPAR(I)=MAX( VEC(3), ERMIN) IC=INDEX(CARD,'/') IF(IC.GT.0) THEN TLBL=CARD(IC+1:IC+10) ELSE TLBL=' ' ENDIF IF(LUBIN.EQ.0) PLBL(I)=TLBL WRITE (LUOUT,480) NPAR,IDPAR(I),PAR(I),ERPAR(I),TLBL 10 NPARTOT=I C C READ IN CORRELATION INFO (OPTIONAL) C 50 NL=NPAR NTOP=NTOP+1-NL*NL IDPAR(0)=NPARTOT WRITE(LUOUT,490) NPARTOT IF(NTOP.LT.0) RETURN IF(NPARTOT.EQ.N) GETPAR=GETVAR(LU,NPAR,VAR(NTOP),LUBIN) IF(GETPAR.LT.0) RETURN IF(GETPAR.EQ.0) THEN C MAKE DEFAULT VAR NL=NTOP DO 60 I=1,NPAR CALL DCOPY(NPAR,ZERO,I0,VAR(NL),I1) 60 NL=NL+NPAR ENDIF NVAL=NPAR NL=NTOP N=NPAR+1 DO 70 I=1,NPARTOT IF(IDPAR(I).GE.0) THEN VAL=ERPAR(I) IF(GETPAR.EQ.0) THEN VAR(NL)=VAL ELSE CALL DSCAL(NVAL,VAL,VAR(NL),NPAR) NVAL=NVAL-1 ENDIF NL=NL+N ENDIF 70 CONTINUE RETURN 100 FORMAT(A) 480 FORMAT(I6,F16.0,E21.13,E15.6,1X,A) 490 FORMAT(I6,' parameters read') END C FUNCTION GETVAR(LU,NPAR,VAR,LUBIN) C C SUBROUTINE TO READ VARIANCE C INTEGER GETVAR,LU,NPAR,LUBIN REAL*8 VAR(NPAR,*) INTEGER I,J,IERR C C READ IN CORRELATION INFO C IF(LUBIN.LE.0) THEN READ (LU,100,IOSTAT=IERR) ((VAR(I,J),J=I,NPAR),I=1,NPAR) IF(IERR.EQ.0) THEN GETVAR=1 RETURN ENDIF ENDIF IF(LUBIN.GT.0) THEN READ (LUBIN,IOSTAT=IERR) ((VAR(I,J),J=I,NPAR),I=1,NPAR) IF(IERR.EQ.0) THEN WRITE(*,*) ' using .unf file for variance' GETVAR=-1 RETURN ENDIF ENDIF GETVAR=0 RETURN 100 FORMAT(8F10.7) END C SUBROUTINE PUTVAR(LU,LUBIN,NPAR,IDPAR,VAR,ERPAR) C C SUBROUTINE TO WRITE VARIANCE C INTEGER LU,LUBIN,NPAR REAL*8 IDPAR(0:*),VAR(NPAR,*),ERPAR(*) REAL*8 VAL INTEGER I,J,K,NDQ,N,NN,IERR NDQ=NPAR C C WRITE CORRELATION INFO C IF(LUBIN.NE.0) THEN WRITE(LUBIN,IOSTAT=IERR) ((VAR(I,J),J=I,NPAR),I=1,NPAR) ENDIF K=0 N=NPAR NN=IDPAR(0) DO 20 I=1,NN IF(IDPAR(I).GE.0) THEN K=K+1 VAL=1./ERPAR(I) CALL DSCAL(N,VAL,VAR(K,K),NDQ) N=N-1 ENDIF 20 CONTINUE WRITE(LU,100) ((VAR(I,J),J=I,NPAR),I=1,NPAR) RETURN 100 FORMAT(8F10.7) END C FUNCTION CALERR(NPAR,VAR,DERV) INTEGER NPAR,K,KK,I1 REAL*8 CALERR,ERR,SUM,DDOT,VAR(0:*),DERV(*) PARAMETER (I1=1) ERR=0.D0 KK=0 DO 1 K=1,NPAR SUM=DDOT(K,DERV,I1,VAR(KK),I1) ERR=ERR+SUM*SUM 1 KK=KK+NPAR CALERR=SQRT(ERR) RETURN END C SUBROUTINE ICOPY(IBGN,IEND,IX,IY) INTEGER IBGN,IEND,I INTEGER*2 IX(*),IY(*) DO 10 I=IBGN,IEND 10 IY(I)=IX(I) RETURN END FUNCTION DEFLIN(IQNFMT,IDQN) INTEGER DEFLIN,IQNFMT,IDQN(12),I,K,KM,NQN NQN=MOD(IQNFMT,10) DO 1 K=1,12 1 IDQN(K)=0 I=MOD(IQNFMT/100,5)+1 KM=1 DO 2 K=I,NQN IDQN(K)=KM IDQN(K+NQN)=IDQN(K)+NQN KM=K 2 CONTINUE DEFLIN=NQN RETURN END FUNCTION GETLIN(LU,IDQN,IQN,XFREQ,XERR,XWT) INTEGER GETLIN,LU,IDQN(12),IQN(12) REAL*8 XFREQ,XERR,XWT C read in .LIN cadr from LU INTEGER J,K,KC,KK,NEG,JD,PCARD,IC,ICH REAL*8 VAL(3) CHARACTER*81 CARD GETLIN=-1 READ (LU,'(A)',IOSTAT=K) CARD IF(K.NE.0) RETURN VAL(1)=0. VAL(2)=0.01 VAL(3)=1. IF(PCARD(CARD(37:81),VAL,3).EQ.0) RETURN XFREQ=VAL(1) XERR=VAL(2) XWT=VAL(3) KC=0 C DECODE QUANTA DO 10 K=1,12 NEG=1 JD=IDQN(K) J=0 DO 11 KK=1,3 KC=KC+1 ICH=ICHAR(CARD(KC:KC)) IC=ICH-ICHAR('0') IF(IC.GE.0.AND.IC.LE.9) THEN J=J*10+IC JD=0 ELSE IF(ICH.EQ.ICHAR('-')) THEN NEG=-1 ENDIF 11 CONTINUE IF(NEG.LT.0) J=-J IF(JD.NE.0) J=IQN(JD) IQN(K)=J IF(J.NE.0) GETLIN=0 10 CONTINUE RETURN END FUNCTION PCARD(CARD,VAL,NVAL) INTEGER PCARD,NVAL CHARACTER*(*) CARD REAL*8 VAL(0:*) C C parses a card for a leading alphabetic character and C NVAL numbers C REAL*8 TMP CHARACTER*1 CH INTEGER*4 ITMP INTEGER MXVAL,NPWR,KVAL,ICH,ICHR,NCHR,NDEC,IPWR,PWRFLG PARAMETER (NPWR=6) LOGICAL NEG,NEWNUM REAL PTEN(0:NPWR) SAVE PTEN DATA PTEN/1.,10.,100.,1000.,10000.,100000.,1000000./ NCHR=LEN(CARD) NEG=.FALSE. NEWNUM=.TRUE. ITMP=0 NDEC=-1 IPWR=0 PWRFLG=0 KVAL=0 MXVAL=NVAL ICHR=0 50 IF(KVAL.LT.MXVAL) THEN ICHR=ICHR+1 IF(ICHR.LE.NCHR) THEN CH=CARD(ICHR:ICHR) ELSE CH=' ' MXVAL=0 ENDIF ICH=ICHAR(CH)-ICHAR('0') IF(ICH.GE.0.AND.ICH.LE.9) THEN C character is a number IF(NDEC.LE.0) THEN ITMP=ICH NDEC=1 ELSE IF(NDEC.EQ.NPWR) THEN C .. now is a good time to convert integer to real IF(NEWNUM) THEN TMP=ITMP NEWNUM=.FALSE. ELSE TMP=TMP*PTEN(NDEC)+ITMP ENDIF ITMP=ICH NDEC=1 ELSE NDEC=NDEC+1 ITMP=ITMP*10+ICH ENDIF IF(PWRFLG.GT.0) IPWR=IPWR-1 ELSE IF(CH.EQ.'.') THEN PWRFLG=1 IF(NDEC.LT.0) NDEC=0 ELSE IF(CH.EQ.'-') THEN NEG=.TRUE. NDEC=0 ELSE IF(CH.EQ.'+') THEN NDEC=0 ELSE C character is not a number or decimal point, +, - IF(NDEC.GE.0) THEN C save results from number decoding IF(PWRFLG.LT.0) THEN C integer follows 'E' IF(NEG) ITMP=-ITMP IPWR=IPWR+ITMP ELSE C finish up number IF(NEWNUM) THEN TMP=ITMP ELSE TMP=ITMP+TMP*PTEN(NDEC) ENDIF IF(NEG) TMP=-TMP ENDIF IF(CH.EQ.'E'.OR.CH.EQ.'e'.OR. : CH.EQ.'D'.OR.CH.EQ.'d') THEN C ..... look for exponent PWRFLG=-1 NDEC=0 ELSE C .. fix up power of 10 in groups of NPWR 10 IF(IPWR.NE.0) THEN IF(IPWR.LT.0) THEN ITMP=-IPWR IF(ITMP.GT.NPWR) ITMP=NPWR TMP=TMP/PTEN(ITMP) IPWR=IPWR+ITMP ELSE ITMP=IPWR IF(ITMP.GT.NPWR) ITMP=NPWR TMP=TMP*PTEN(ITMP) IPWR=IPWR-ITMP IF(ABS(TMP).GT.1.E+32) IPWR=0 ENDIF GO TO 10 ENDIF VAL(KVAL)=TMP KVAL=KVAL+1 PWRFLG=0 NDEC=-1 IPWR=0 ENDIF NEG=.FALSE. NEWNUM=.TRUE. ITMP=0 ENDIF IF(CH.EQ.',') THEN IF(NDEC.EQ.-2) KVAL=KVAL+1 NDEC=-2 ELSE IF(CH.EQ.'/') THEN C end of line character KVAL=MXVAL ENDIF ENDIF GO TO 50 ENDIF PCARD=KVAL RETURN END