C--------------------------------------------------- C Spectrum of perturbations: P(k), sigma_8,V(50h-1Mpc) C Anatoly Klypin December 1997 C Holtzman appoximation C or BBKS-style + Hu&Sugiyama approx C Needs cdm.fit file -- parameters of models and fitting C----------------------------------------------------- C Om = Omega_0 Omc = Omega_cdm_present C Omb = Omega_baryon Omnu =Omega_neutrino_present C Ocurv = Omega_curvature_present C hsmall = H_0/100km/s/Mpc = the Hubble constant C ns = slope of the power spectum (ns=1 for Harr-Zeld) C AEXPN = expansion parameter = 1/(1+z), z =redshift C All calculations are done in real Mpc, C Only final outputs are scaled to h^{-1} C----------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NtabM = 100000) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall COMMON / Coeff/ deltac22,deltac,rf,Uklow,Ukup COMMON /PTABLE/xkt(0:NtabM),Pkt(0:NtabM),StepK,alog0,iFlag,Ntab common /sigmx/sig,ai CHARACTER Name*16,Answer*1,Header*45,StartFile*16 REAL*8 INTG,ns,INPUT,P2,P,Ptophat,Pgauss,Vtophat EXTERNAL INTG,P2,P,Ptophat,Pgauss,Vtophat DATA StartFile/'InStart.dat'/ WRITE (*,*) '================< Enter File Name for output file: ' READ (*,'(a)') Name OPEN(1,file=Name) OPEN(2,file='cdm.fit') ! parameters of models C Constants PI =3.1415926535 AEXPN =1. WRITE (*,'(A,$)') ' Enter sigma_8, slope n, SigmaX: ' READ (*,*) Sigma8,ns,sig ai = 1.d-3 CALL MODEL(Line) H =100.*hsmall ! Hubble constant DO a=0.1,1.,0.1 CALL AGE(t0,GrowthDen,GrowthVel,a) ENDDO AEXPN =1. a =AEXPN write (*,*) ' a=',a,' Aexpn=',AEXPN CALL AGE(t0,GrowthDen_0,GrowthVel_0,a) c Do i=1,Ntab,100 | test interpolation c xx = xkt(i)*1.12 c yy = Ppk(xx) c EndDo c Stop rf =8./hsmall ! top-hat for sigma_8 Sn = (Sigma8)**2 / & INTG(Ptophat,1.d-5,2.d+1) ! normalization of P(k) c bias parameter Sigma8 = sqrt(Sn*( & INTG(Ptophat,1.d-5,2.d+1))) write (*,10) hsmall,Sigma8,ns ! check if all is self-consistent write (1,10) hsmall,Sigma8,ns 10 Format(' Hubble=',f7.3,' Sigma_8 =',G11.3, & ' Slope n=',f7.2) rf =50./hsmall ! top-hat for bulk velocity Veloc = GrowthVel*H* & sqrt(Sn*(INTG(Vtophat,1.d-5,1.d+0))) write (*,30) Veloc,rf*hsmall write (1,30) Veloc,rf*hsmall 30 Format(' Bulk Velocity =',F8.2, . 'km/s for radius top-hat=',f8.1,'h-1Mpc') write (1,*) ' k/h Pk*h^3 Power Spectrum at z=',1/AEXPN-1. DO i=1,110 w =1.e-4*10.**(i/20.) Pk0 =P(w)*Sn*hsmall**3*(2.*pi**2) write (1,14) w/hsmall,Pk0 ENDDO 14 Format(2G11.4,5G12.4) C..........................................................Power spectrum and C Amplitude of fluctuations in a Box at redshift z WRITE (*,'(A,$)') & 'Enter Box(h^-1Mpc), Nparticles(1D) and redshift:' READ (*,*) Box, NROW, z Box =Box/hsmall ! scale it to real Mpc Uklow =2.*pi/Box ! frequency range for integrals Ukup =2.*pi/Box*(NROW/2.) AEXPN =1./(1.+z) ! expansion parameter a =AEXPN CALL AGE(t0,GrowthDen,GrowthVel,a) sigma =(GrowthDen/GrowthDen_0)* & sqrt(Sn*INTG(P2,Uklow/sqrt(2.),Ukup)) WRITE (*,40) z,sigma,NROW,Box*hsmall 40 format(' z=',f8.3,' delta\rho/rho in box=',f9.5,/ . 5X,'Particles=',i4,' Box=',f8.2) write (1,*) ' k/h Pk*h^3 Power Spectrum at z=',1/AEXPN-1. DO i=1,130 w =1.e-4*10.**(i/20.) Pk0 =(GrowthDen/GrowthDen_0)**2*P(w)*Sn*hsmall**3*(2.*pi**2) write (1,14) w/hsmall,Pk0 ENDDO C............................................................................................ WRITE (*,'(A,$)') ' Do you need a file to run PMstart (Yes/No) ' READ (*,'(a)') Answer IF(Answer.eq.'Y'.or.Answer.eq.'y')Then ! prepare input file OPEN(30,file=StartFile) ! parameters of a run write(30,'(A)') 'Yes Read parameters from cdm.fit' HEADER='N=128x256 L=20h-1CDMsig=0.239 z=30-----------' write (*,*) '------ Enter Header up to 45 characters' write (*,*) ' Example:' write (*,'(A)') HEADER read (*,'(A)') HEADER write (30,'(A)') HEADER write (30,*) AEXPN,' Expansion Parameter' ASTEP =INPUT(' Step of the expansion parameter =') write (30,*) ASTEP,' Step in dAEXPN ' write (30,*) sigma,' DRho/rho in box ' write (30,*) ns, ' Slope of P(k) ' write (30,*) Box*hsmall,' Box in Mpc/h ' write (30,*) Ocurv,' Omega_curvature ' Nseed =INPUT(' Enter random seed integer 1-2^31 =') write (30,*) Nseed,' Random seed ' write (30,*) Line,' Line number in cdm.fit' Write (*,'(A,$)') ' Enter Min and Max number of mass species =' read (*,*) Nmin_lev,Nmax_lev write (30,*) Nmin_lev,Nmax_lev,' Min and Max number of species' write (*,*) ' Results were written to file: ',StartFile CLOSE (30) ENDIF STOP END c----------------------------------- Pcold(k)*k^2 REAL*8 FUNCTION P2(WK) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall COMMON / Coeff/ deltac22,deltac,rf,Uklow,Ukup Real*8 ns P2=WK**2*P(WK) RETURN END C--------------------------------------------- C Power spectrum C wk = k= in real Mpc REAL*8 FUNCTION P(wk) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NtabM = 100000) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall COMMON / FACT/ qqscaleb COMMON /PTABLE/xkt(0:NtabM),Pkt(0:NtabM),StepK,alog0,iFlag,Ntab Real*8 ns If(iFlag.gt.0)Then ! use approximations If(Par(6).ne.0.)Then ! Holtzman approx sk= sqrt(wk) P = wk**ns / . (1.+sk*(Par(2) . +sk*(Par(3) . +sk*(Par(4) . +sk*(Par(5) )))) )**(2.*Par(6)) Else ! BBKS + Sugiyama approx c Gammaeff =Om*hsmall/exp(Omb*(1.+sqrt(hsmall/0.5)/Om)) c Q = wk/hsmall/Gammaeff Q = wk*qqscaleb P = wk**ns*( LOG(1.+Par(1)*Q)/(Par(1)*Q) )**2/ . sqrt(1.+Q*(Par(2)+ . Q*(Par(3)**2+ . Q*(Par(4)**3+ . Q*(Par(5)**4) )))) EndIf Else ! use table for P(k) P = Ppk(wk) EndIf RETURN END C------------------------------------------------- C Age of the Universe: t0 (z=0) C Expansion parameter: a =1/(1+z) C Growth_Rate_Density at a: GrowthDen C GrowthDen =\delta\rho/\rho C normalized to GrowthDen =a for Omega=1 C GrowthDen < 1 at a=1 for Lcdm or OCDM C Growth_Rate_velocity at a: GrowthVel C GrowthVel = V_peculiar= (a/delta)(d delta/d a) SUBROUTINE AGE(t0,GrowthDen,GrowthVel,a) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (zero =1.d-12) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall Common /OmegRat/ OmLOm0,OmcOm0 common /sigmx/sig,ai Real*8 INTG,ns External Hnorm,Hage,Hgrow Oml =Max(Abs(1.-Om-Ocurv),zero) OmLOm0 =Oml/Om OmcOm0 =Ocurv/Om t0 =9.766/hsmall/sqrt(Om)*INTG(Hage,zero,a) z =1./a-1. Hubble = hsmall*sqrt(Om/a**3)*Hnorm(a) c ww =INTG(Hgrow,zero,ai) wwa =INTG(Hgrow,zero,a) GrowthVel = -(1.5+OmcOm0*a)/Hnorm(a)**2+ & sqrt(a**5)/Hnorm(a)**3/wwa GrowthDeni =2.5*Hnorm(ai)/sqrt(ai**3)*ww GrowthDen=grdc(a,oml,om,sig,ai,growthdeni) c write (*,20) t0,z,a,Hubble,GrowthDen,GrowthVel write (1,20) t0,z,a,Hubble,GrowthDen,GrowthVel 20 format(' Age=',f8.3,' z=',f6.2,' a=',f6.3, . ' Hubble/100=',f8.3,' GrowthRateDen=',f7.4, . ' GrowthRateVelocity=',f7.4) Return End c----------------------------------- P(k)*k^2*gauss(rf) REAL*8 FUNCTION Pgauss(WK) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall COMMON / Coeff/ deltac22,deltac,rf,Uklow,Ukup Real*8 ns Pgauss=WK**2*P(WK)*EXP(-(rf*WK)**2) RETURN END c----------------------------------- P(k)*k^4*gauss(rf) REAL*8 FUNCTION P4gauss(WK) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall COMMON / Coeff/ deltac22,deltac,rf,Uklow,Ukup Real*8 ns P4gauss=WK**4*P(WK)*EXP(-(rf*WK)**2) RETURN END C-------------------------------------- P*k^2*Top-Hat Filter REAL*8 FUNCTION Ptophat(wk) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall COMMON / Coeff/ deltac22,deltac,rf,Uklow,Ukup Real*8 ns IF (wk.lt.1.d-4) THEN Ptophat =wk**ns*wk**2 ELSE X =wk*rf TOPHAT =( (SIN(X)-x*COS(X))*3./X**3 )**2 Ptophat=P(wk)*wk**2*TOPHAT ENDIF RETURN END C-------------------------------------- P*Top-Hat Filter REAL*8 FUNCTION Vtophat(wk) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall COMMON / Coeff/ deltac22,deltac,rf,Uklow,Ukup Real*8 ns IF (wk.lt.1.d-4) THEN Vtophat =wk**ns ELSE X=wk*rf TOPHAT=( (SIN(X)-x*COS(X))*3./X**3 )**2 Vtophat=P(wk)*TOPHAT ENDIF RETURN END C-------------------------------------- Linear Pair-wise V REAL*8 FUNCTION Vpairwise(wk) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall COMMON / Coeff/ deltac22,deltac,rf,Uklow,Ukup Real*8 ns X =wk*rf Vpairwise=P(wk)*(1.-sin(X)/X) RETURN END C________________________________________Read parameters of the model C Om = Omega_0 C Omb = Omega_baryon C Omnu = Omega_neutrino C hsmall = hubble constant C Par = fitting parameters C Par(6) = 0 --- bbks-style+Hu&Sugiyama C Par(6) ne 0 --- holtzman-style SUBROUTINE MODEL(Line) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NtabM = 100000) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall COMMON / FACT/ qqscaleb COMMON /PTABLE/xkt(0:NtabM),Pkt(0:NtabM),StepK,alog0,iFlag,Ntab Real*8 ns,INPUT Character Header*79 Line =INPUT(' Enter Line Number in cdm.fit= ') If(Line.gt.0)Then iFlag = Line Read(2,'(A)') Header If(Line.gt.2)Then Do i=1,Line -2 Read (2,*) a EndDo EndIf Read (2,*) Om,Omb,Omc,Omnu,hsmall,Par Ocurv =INPUT(' Enter Omega curvature at z=0 =>') theta = 2.726/2.7 ! = T_cmb/2.7 Ob0 =Omb/Om ! Hu&Sugiyama fitting Omh2 =Om*hsmall**2 a1 =(46.9*Omh2)**0.670*(1.+(32.1*Omh2)**(-0.532)) a2 =(12.0*Omh2)**0.424*(1.+(45.*Omh2)**(-0.582)) alpha =1./a1**Ob0/a2**(Ob0**3) qqscaleb = theta**2/(1.-Ob0)**0.60/sqrt(alpha)/Omh2 Write (*,20)Om,Omb,Omc,Omnu,hsmall,Par Write (1,20)Om,Omb,Omc,Omnu,hsmall,Par 20 Format(' Model: Om_0=',F5.3,' Om_baryon=',F5.3, . ' Om_cold=',F5.3,' O_nu=',F5.3, . ' hsmall=',F4.2,/8x,'Parameters=',(6G9.3)) Else iFlag = -1 hsmall = 0.7 Omb = 0.024/hsmall**2 Om = 0.3 Ocurv = 0. Omc = Om CALL ReadTable Write (*,22)Om,Omb,hsmall Write (1,22)Om,Omb,hsmall 22 Format(' Model: Om_0=',F5.3,' Om_baryon=',F5.3, . ' hsmall=',F4.2, ' <== Using table for P(k)') EndIf Return End C-------------------------------------------- Simpson integration REAL*8 FUNCTION INTG(FUNC,A,B) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (EPS=5.0d-4, JMAX=22) EXTERNAL FUNC OST=-1.D30 OS= -1.D30 ST =0. DO 11 J=1,JMAX CALL TRAPZD(FUNC,A,B,ST,J) INTG=(4.0d0*ST-OST)/3.0d0 IF (ABS(INTG-OS).Le.EPS*ABS(OS)) RETURN OS=INTG OST=ST 11 CONTINUE WRITE (*,*)'Integration did not converge' RETURN END C---------------------------------------------- SUBROUTINE TRAPZD(FUNCC,A,B,S,N) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) SAVE IT EXTERNAL FUNCC IF (N.EQ.1) THEN S=0.5d0*(B-A)*(FUNCC(A)+FUNCC(B)) IT=1 ELSE TNM=IT DEL=(B-A)/TNM X=A+0.5D0*DEL SUM=0.0D0 DO 11 J=1,IT SUM=SUM+FUNCC(X) X=X+DEL 11 CONTINUE S=0.5D0*(S+(B-A)*SUM/TNM) IT=2*IT ENDIF RETURN END C-------------------------------------------- Simpson integration REAL*8 FUNCTION INTG1(FUNC,A,B) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (EPS=1.0d-3, JMAX=22) EXTERNAL FUNC OST=-1.d30 OS= -1.d30 ST =0. DO 11 J=1,JMAX CALL TRAPZD1(FUNC,A,B,ST,J) INTG1=(4.0D0*ST-OST)/3.0D0 IF (ABS(INTG1-OS).Le.EPS*ABS(OS)) RETURN OS=INTG1 OST=ST 11 CONTINUE c WRITE (1,*)'Integration did not converge' RETURN END C---------------------------------------------- SUBROUTINE TRAPZD1(FUNCC,A,B,S,N) C--------------------------------------- IMPLICIT Real*8 (A-H,O-Z) SAVE IT EXTERNAL FUNCC IF (N.EQ.1) THEN S=0.5D0*(B-A)*(FUNCC(A)+FUNCC(B)) IT=1 ELSE TNM=IT DEL=(B-A)/TNM X=A+0.5D0*DEL SUM=0.0D0 DO 11 J=1,IT SUM=SUM+FUNCC(X) X=X+DEL 11 CONTINUE S=0.5D0*(S+(B-A)*SUM/TNM) IT=2*IT ENDIF RETURN END C---------------------------------- Read in variables REAL*8 FUNCTION INPUT(text) C----------------------------------------- Character text*(*) write (*,'(A,$)')text read (*,*) x INPUT =x Return End C----------------------------------------- REAL*8 FUNCTION sinhh(x) C----------------------------------------- IMPLICIT REAL*8 (A-H,P-Z) sinhh =0.5*(exp(MIN(x,33.d+0))-exp(-MIN(x,33.d+0))) RETURN END C--------------------------------------- REAL*8 FUNCTION Hh(x) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) Hh =(sqrt(x/(1.d+0 +x**3)))**3 Return End C--------------------------------------- REAL*8 FUNCTION Hnorm(x) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) Common /OmegRat/ OmLOm0,OmcOm0 Hnorm =sqrt(1.d+0 & +OmLOm0*x**3 & +OmcOm0*x ) Return End C--------------------------------------- REAL*8 FUNCTION Hage(x) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) Hage =sqrt(x)/Hnorm(x) Return End C--------------------------------------- REAL*8 FUNCTION Hgrow(x) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) Hgrow =(sqrt(x)/Hnorm(x))**3 Return End C--------------------------------------- REAL*8 FUNCTION Ppk(x) c interpolate table with p(k) c x is in real 1/Mpc C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NtabM = 100000) COMMON /PTABLE/xkt(0:NtabM),Pkt(0:NtabM),StepK,alog0,iFlag,Ntab If(x.ge.xkt(Ntab))Then ! slope is ns =-3 Ppk =Pkt(Ntab)/(x/xkt(Ntab))**3 Return EndIf If(x.lt.xkt(1))Then ! slope is ns=1 Ppk =Pkt(1)*(x/xkt(1)) Return EndIf ind = INT((log10(x)-alog0)/StepK) +1 dk = xkt(ind+1)-xkt(ind) Ppk = (Pkt(ind)*(xkt(ind+1)-x)+Pkt(ind+1)*(x-xkt(ind)))/dk c write (*,'(5x," k=",g12.4," table=",2g12.4," P=",3g12.4)') c & x,xkt(ind),xkt(ind+1), c & Ppk,Pkt(ind),Pkt(ind+1) Return End C--------------------------------------- SUBROUTINE ReadTable c interpolate table with p(k) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NtabM = 100000) COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Ocurv,Par(6),AEXPN,Sn,ns,hsmall COMMON /PTABLE/xkt(0:NtabM),Pkt(0:NtabM),StepK,alog0,iFlag,Ntab REAL*8 ns Open(50,file='WHUpk.dat') twop = 2.*3.14145926**2 Ntab =0 10 read(50,*,end=30,err=30)xx,pp Ntab =Ntab +1 xkt(Ntab) =xx *hsmall ! scale to real Mpc Pkt(Ntab) = pp /xx**3*twop goto 10 30 write (*,*) ' Read ',Ntab,' lines from file WHUpk.dat' If(Ntab.le.1)stop 'wrong table for p(k)' StepK =log10(xkt(2)/xkt(1)) alog0 = log10(xkt(1)) Do k =2,Ntab-1 ! test that spacing is the same ss =log10(xkt(k+1)/xkt(k)) If(abs(ss/StepK-1.).gt.1.e-3)Then write (*,*) ' error in K spacing. k=',k,xkt(k+1),xkt(k) STOP EndIf EndDo Return End c function grdc(a,omlax,omemx,sigx,ai,dci) c ********************************** c growth of density contrast c ******************************* implicit real*8 (a-h,o-z) common /par1/omem,omla,sig common /par2/j0 dimension yy(2),info(15),rwork(162),iwork(33) dimension rpar(1) data lrw,liw/162,33/ external df,podint c j0=1 omla=omlax omem=omemx sig=sigx c c ********* pocetni uvjeti *************** c derdci=dci/ai c c ********** integracija ******************* c info(1)=0 info(2)=0 info(3)=0 info(4)=0 rtol=1.d-7 atol=1.d-7 xx=ai xxout=a yy(1)=dci yy(2)=derdci 2 call ddeabm(df,2,xx,yy,xxout,info,rtol,atol,idid,rwork, 1lrw,iwork,liw,rpar,ipar) if(idid.eq.-2)then write(6,*)'*****************' write(6,*)'message from grdc and solver:' write(6,*)'changed tolerance' write(6,*)'rtol=',rtol,' atol=',atol write(6,*)'xx=',xx,' xxout=',xxout write(6,*)'******************' info(1)=1 goto 2 endif if(idid.lt.0)then write(6,*)'****************' write(6,*)'message from grdc and solver:' write(6,*)'idid=',idid write(6,*)'xx=',xx,' xxout=',xxout write(6,*)'*****************' endif grdc=yy(1) return end c c c-------------------------------------------------------------------- subroutine df(x,u,uprime,rpar,ipar) implicit real*8 (a-h,o-z) common /par1/omem,omla,sig common /par2/j0 dimension u(1),uprime(1),rpar(1) c pom=omem/x**3+omla c uprime(1)=u(2) if(j0.eq.2)goto 1 uprime(2)=-1.5d0/x*(omem/x**3+2.d0*omla)/pom*u(2) 1+1.5d0/x**5*omem/pom*u(1)+1.5d0*omem**(4.d0/3.d0) 2*dsqrt(sig)*x**(-6.d0-0.25d0)*pom**(-1.5d0)*u(1) goto 2 1 continue uprime(2)=-1.5d0/x*(omem/x**3+2.d0*omla)/pom*u(2) 1+1.5d0/x**5*omem/pom*u(1)+1.5d0*omem**(4.d0/3.d0) 2*dsqrt(sig)*x**(-6.d0-0.25d0)*pom**(-1.5d0)*u(1) 3+0.25d0*dsqrt(sig)*omem**(1.d0/3.d0)*x**(-3.d0-0.25d0) 4/dsqrt(pom)*u(1) 2 continue return end c c c-------------------------------------------------------------------- function bd(a) implicit real*8 (a-h,o-z) common /par1/omem,omla,sig common /par2/j0 if(j0.eq.1)goto 1 if(j0.eq.2)goto 2 1 bd=dsqrt(sig)*omem**(-1.d0/6.d0)*dsqrt(a) goto 3 2 bd=7.d0/6.d0*dsqrt(sig)*omem**(-1.d0/6.d0)*a**(0.25d0) 3 continue return end c-------------------------------------------------------------------- c function alfa(a) implicit real*8 (a-h,o-z) alfa=-1.d0/4.d0+0.5d0*dsqrt(0.25d0+6.d0*(1.d0+bd(a))) return end c c-------------------------------------------------------------------- c c function h(a) c implicit real*8 (a-h,o-z) c common /par1/omem,omla,sig c common /par2/j0 c h=dsqrt(omem/a**3+omla) c return c end c c c-------------------------------------------------------------------- function podint(x) implicit real*8 (a-h,o-z) podint=x**(1.5d0)/(1.d0+x**3)**1.5d0 return end c-------------------------------------------------------------------- SUBROUTINE DDEABM(DF,NEQ,T,Y,TOUT,INFO,RTOL,ATOL,IDID,RWORK,LRW, 1 IWORK,LIW,RPAR,IPAR) C***BEGIN PROLOGUE DDEABM C***DATE WRITTEN 820301 (YYMMDD) C***REVISION DATE 830701 (YYMMDD) C***CATEGORY NO. I1A1B C***KEYWORDS ADAMS METHOD,DEPAC,DOUBLE PRECISION, C INITIAL VALUE PROBLEMS,ODE, C ORDINARY DIFFERENTIAL EQUATIONS,PREDICTOR-CORRECTOR C***AUTHOR SHAMPINE, L.F. C WATTS, H.A. C***PURPOSE Solves initial value problems in ordinary differential C equations using an Adams-Bashforth method. C***DESCRIPTION C C This is the Adams code in the package of differential equation C solvers DEPAC, consisting of the codes DDERKF, DDEABM, and DDEBDF. C Design of the package was by L. F. Shampine and H. A. Watts. C It is documented in C SAND79-2374 , DEPAC - Design of a User Oriented Package of ODE C Solvers. C DDEABM is a driver for a modification of the code ODE written by C L. F. Shampine and M. K. Gordon C Sandia Laboratories C Albuquerque, New Mexico 87185 C C ********************************************************************** C * ABSTRACT * C ************ C C Subroutine DDEABM uses the Adams-Bashforth-Moulton C Predictor-Corrector formulas of orders one through twelve to C integrate a system of NEQ first order ordinary differential C equations of the form C DU/DX = DF(X,U) C when the vector Y(*) of initial values for U(*) at X=T is given. C The subroutine integrates from T to TOUT. It is easy to continue the C integration to get results at additional TOUT. This is the interval C mode of operation. It is also easy for the routine to return with C the solution at each intermediate step on the way to TOUT. This is C the intermediate-output mode of operation. C C DDEABM uses subprograms DDE2, DSTEP2, DINTRP, DHSTRT, DVNORM, C D1MACH, and the error handling routine XERRWV. The only machine C dependent parameters to be assigned appear in D1MACH. C C ********************************************************************** C * Description of The Arguments To DDEABM (An Overview) * C ********************************************************************** C C The Parameters are C C DF -- This is the name of a subroutine which you provide to C define the differential equations. C C NEQ -- This is the number of (first order) differential C equations to be integrated. C C T -- This is a value of the independent variable. C C Y(*) -- This array contains the solution components at T. C C TOUT -- This is a point at which a solution is desired. C C INFO(*) -- The basic task of the code is to integrate the C differential equations from T to TOUT and return an C answer at TOUT. INFO(*) is an integer array which is used C to communicate exactly how you want this task to be C carried out. C C RTOL, ATOL -- These quantities represent relative and absolute C error tolerances which you provide to indicate how C accurately you wish the solution to be computed. You may C choose them to be both scalars or else both vectors. C C IDID -- This scalar quantity is an indicator reporting what C the code did. You must monitor this integer variable to C decide what action to take next. C C RWORK(*), LRW -- RWORK(*) is a DOUBLE PRECISION work array of C length LRW which provides the code with needed storage C space. C C IWORK(*), LIW -- IWORK(*) is an INTEGER work array of length LIW C which provides the code with needed storage space. C C RPAR, IPAR -- These are DOUBLE PRECISION and INTEGER parameter C arrays which you can use for communication between your C calling program and the DF subroutine. C C Quantities which are used as input items are C NEQ, T, Y(*), TOUT, INFO(*), C RTOL, ATOL, RWORK(1), LRW and LIW. C C Quantities which may be altered by the code are C T, Y(*), INFO(1), RTOL, ATOL, C IDID, RWORK(*) and IWORK(*). C C ********************************************************************** C * INPUT -- What To Do On The First Call To DDEABM * C ********************************************************************** C C The first call of the code is defined to be the start of each new C problem. Read through the descriptions of all the following items, C provide sufficient storage space for designated arrays, set C appropriate variables for the initialization of the problem, and C give information about how you want the problem to be solved. C C C DF -- Provide a subroutine of the form C DF(X,U,UPRIME,RPAR,IPAR) C to define the system of first order differential equations C which is to be solved. For the given values of X and the C vector U(*)=(U(1),U(2),...,U(NEQ)) , the subroutine must C evaluate the NEQ components of the system of differential C equations DU/DX=DF(X,U) and store the derivatives in the C array UPRIME(*), that is, UPRIME(I) = * DU(I)/DX * for C equations I=1,...,NEQ. C C Subroutine DF must NOT alter X or U(*). You must declare C the name df in an external statement in your program that C calls DDEABM. You must dimension U and UPRIME in DF. C C RPAR and IPAR are DOUBLE PRECISION and INTEGER parameter C arrays which you can use for communication between your C calling program and subroutine DF. They are not used or C altered by DDEABM. If you do not need RPAR or IPAR, C ignore these parameters by treating them as dummy C arguments. If you do choose to use them, dimension them in C your calling program and in DF as arrays of appropriate C length. C C NEQ -- Set it to the number of differential equations. C (NEQ .GE. 1) C C T -- Set it to the initial point of the integration. C You must use a program variable for T because the code C changes its value. C C Y(*) -- Set this vector to the initial values of the NEQ solution C components at the initial point. You must dimension Y at C least NEQ in your calling program. C C TOUT -- Set it to the first point at which a solution C is desired. You can take TOUT = T, in which case the code C will evaluate the derivative of the solution at T and C return. Integration either forward in T (TOUT .GT. T) or C backward in T (TOUT .LT. T) is permitted. C C The code advances the solution from T to Tout using C step sizes which are automatically selected so as to C achieve the desired accuracy. If you wish, the code will C return with the solution and its derivative following C each intermediate step (intermediate-output mode) so that C you can monitor them, but you still must provide TOUT in C accord with the basic aim of the code. C C The first step taken by the code is a critical one C because it must reflect how fast the solution changes near C the initial point. The code automatically selects an C initial step size which is practically always suitable for C the problem. By using the fact that the code will not step C past TOUT in the first step, you could, if necessary, C restrict the length of the initial step size. C C For some problems it may not be permissible to integrate C past a point TSTOP because a discontinuity occurs there C or the solution or its derivative is not defined beyond C TSTOP. When you have declared a TSTOP point (see INFO(4) C and RWORK(1)), you have told the code not to integrate C past TSTOP. In this case any TOUT beyond TSTOP is invalid C input. C C INFO(*) -- Use the INFO array to give the code more details about C how you want your problem solved. This array should be C dimensioned of length 15 to accomodate other members of C DEPAC or possible future extensions, though DDEABM uses C only the first four entries. You must respond to all of C the following items which are arranged as questions. The C simplest use of the code corresponds to answering all C questions as YES ,i.e. setting ALL entries of INFO to 0. C C INFO(1) -- This parameter enables the code to initialize C itself. You must set it to indicate the start of every C new problem. C C **** Is this the first call for this problem ... C YES -- set INFO(1) = 0 C NO -- not applicable here. C See below for continuation calls. **** C C INFO(2) -- How much accuracy you want of your solution C is specified by the error tolerances RTOL and ATOL. C The simplest use is to take them both to be scalars. C To obtain more flexibility, they can both be vectors. C The code must be told your choice. C C **** Are both error tolerances RTOL, ATOL scalars ... C YES -- set INFO(2) = 0 C and input scalars for both RTOL and ATOL C NO -- set INFO(2) = 1 C and input arrays for both RTOL and ATOL **** C C INFO(3) -- The code integrates from T in the direction C of TOUT by steps. If you wish, it will return the C computed solution and derivative at the next C intermediate step (the intermediate-output mode) or C TOUT, whichever comes first. This is a good way to C proceed if you want to see the behavior of the solution. C If you must have solutions at a great many specific C TOUT points, this code will compute them efficiently. C C **** Do you want the solution only at C TOUT (and not at the next intermediate step) ... C YES -- set INFO(3) = 0 C NO -- set INFO(3) = 1 **** C C INFO(4) -- To handle solutions at a great many specific C values TOUT efficiently, this code may integrate past C TOUT and interpolate to obtain the result at TOUT. C Sometimes it is not possible to integrate beyond some C point TSTOP because the equation changes there or it is C not defined past TSTOP. Then you must tell the code C not to go past. C C **** Can the integration be carried out without any C Restrictions on the independent variable T ... C YES -- set INFO(4)=0 C NO -- set INFO(4)=1 C and define the stopping point TSTOP by C setting RWORK(1)=TSTOP **** C C RTOL, ATOL -- You must assign relative (RTOL) and absolute (ATOL) C error tolerances to tell the code how accurately you want C the solution to be computed. They must be defined as C program variables because the code may change them. You C have two choices -- C Both RTOL and ATOL are scalars. (INFO(2)=0) C Both RTOL and ATOL are vectors. (INFO(2)=1) C In either case all components must be non-negative. C C The tolerances are used by the code in a local error test C at each step which requires roughly that C ABS(LOCAL ERROR) .LE. RTOL*ABS(Y)+ATOL C for each vector component. C (More specifically, a Euclidean norm is used to measure C the size of vectors, and the error test uses the magnitude C of the solution at the beginning of the step.) C C The true (global) error is the difference between the true C solution of the initial value problem and the computed C approximation. Practically all present day codes, C including this one, control the local error at each step C and do not even attempt to control the global error C directly. Roughly speaking, they produce a solution Y(T) C which satisfies the differential equations with a C residual R(T), DY(T)/DT = DF(T,Y(T)) + R(T) , C and, almost always, R(T) is bounded by the error C tolerances. Usually, but not always, the true accuracy of C the computed Y is comparable to the error tolerances. This C code will usually, but not always, deliver a more accurate C solution if you reduce the tolerances and integrate again. C By comparing two such solutions you can get a fairly C reliable idea of the true error in the solution at the C bigger tolerances. C C Setting ATOL=0. results in a pure relative error test on C that component. Setting RTOL=0. results in a pure absolute C error test on that component. A mixed test with non-zero C RTOL and ATOL corresponds roughly to a relative error C test when the solution component is much bigger than ATOL C and to an absolute error test when the solution component C is smaller than the threshold ATOL. C C Proper selection of the absolute error control parameters C ATOL requires you to have some idea of the scale of the C solution components. To acquire this information may mean C that you will have to solve the problem more than once. In C the absence of scale information, you should ask for some C relative accuracy in all the components (by setting RTOL C values non-zero) and perhaps impose extremely small C absolute error tolerances to protect against the danger of C a solution component becoming zero. C C The code will not attempt to compute a solution at an C accuracy unreasonable for the machine being used. It will C advise you if you ask for too much accuracy and inform C you as to the maximum accuracy it believes possible. C C RWORK(*) -- Dimension this DOUBLE PRECISION work array of length C LRW in your calling program. C C RWORK(1) -- If you have set INFO(4)=0, you can ignore this C optional input parameter. Otherwise you must define a C stopping point TSTOP by setting RWORK(1) = TSTOP. C (For some problems it may not be permissible to integrate C past a point TSTOP because a discontinuity occurs there C or the solution or its derivative is not defined beyond C TSTOP.) C C LRW -- Set it to the declared length of the RWORK array. C You must have LRW .GE. 120+21*NEQ C C IWORK(*) -- Dimension this INTEGER work array of length LIW in C your calling program. C C LIW -- Set it to the declared length of the IWORK array. C You must have LIW .GE. 33 C C RPAR, IPAR -- These are parameter arrays, of DOUBLE PRECISION and C INTEGER type, respectively. you can use them for C communication between your program that calls DDEABM and C the DF subroutine. They are not used or altered by C DDEABM. If you do not need RPAR or IPAR, ignore these C parameters by treating them as dummy arguments. If you do C choose to use them, dimension them in your calling program C and in DF as arrays of appropriate length. C C ********************************************************************** C * OUTPUT -- After Any Return From DDEABM * C ********************************************************************** C C The principal aim of the code is to return a computed solution at C TOUT, although it is also possible to obtain intermediate results C along the way. To find out whether the code achieved its goal C or if the integration process was interrupted before the task was C completed, you must check the IDID parameter. C C C T -- The solution was successfully advanced to the C output value of T. C C Y(*) -- Contains the computed solution approximation at T. C You may also be interested in the approximate derivative C of the solution at T. It is contained in C RWORK(21),...,RWORK(20+NEQ). C C IDID -- Reports what the code did C C *** Task Completed *** C Reported by positive values of IDID C C IDID = 1 -- A step was successfully taken in the C intermediate-output mode. The code has not C yet reached TOUT. C C IDID = 2 -- The integration to TOUT was successfully C completed (T=TOUT) by stepping exactly to TOUT. C C IDID = 3 -- The integration to TOUT was successfully C completed (T=TOUT) by stepping past TOUT. C Y(*) is obtained by interpolation. C C *** Task Interrupted *** C Reported by negative values of IDID C C IDID = -1 -- A large amount of work has been expended. C (500 steps attempted) C C IDID = -2 -- The error tolerances are too stringent. C C IDID = -3 -- The local error test cannot be satisfied C because you specified a zero component in ATOL C and the corresponding computed solution C component is zero. Thus, a pure relative error C test is impossible for this component. C C IDID = -4 -- The problem appears to be stiff. C C IDID = -5,-6,-7,..,-32 -- Not applicable for this code C but used by other members of DEPAC or possible C future extensions. C C *** Task Terminated *** C Reported by the value of IDID=-33 C C IDID = -33 -- The code has encountered trouble from which C it cannot recover. A message is printed C explaining the trouble and control is returned C to the calling program. For example, this occurs C when invalid input is detected. C C RTOL, ATOL -- These quantities remain unchanged except when C IDID = -2. In this case, the error tolerances have been C increased by the code to values which are estimated to be C appropriate for continuing the integration. However, the C reported solution at T was obtained using the input values C of RTOL and ATOL. C C RWORK, IWORK -- Contain information which is usually of no C interest to the user but necessary for subsequent calls. C However, you may find use for C C RWORK(11)--which contains the step size H to be C attempted on the next step. C C RWORK(12)--if the tolerances have been increased by the C code (IDID = -2) , they were multiplied by the C value in RWORK(12). C C RWORK(13)--which contains the current value of the C independent variable, i.e. the farthest point C integration has reached. this will be different C from T only when interpolation has been C performed (IDID=3). C C RWORK(20+I)--which contains the approximate derivative C of the solution component Y(I). In DDEABM, it C is obtained by calling subroutine DF to C evaluate the differential equation using T and C Y(*) when IDID=1 or 2, and by interpolation C when IDID=3. C C ********************************************************************** C * INPUT -- What To Do To Continue The Integration * C * (calls after the first) * C ********************************************************************** C C This code is organized so that subsequent calls to continue the C integration involve little (if any) additional effort on your C part. You must monitor the IDID parameter in order to determine C what to do next. C C Recalling that the principal task of the code is to integrate C from T to TOUT (the interval mode), usually all you will need C to do is specify a new TOUT upon reaching the current TOUT. C C Do not alter any quantity not specifically permitted below, C in particular do not alter NEQ, T, Y(*), RWORK(*), IWORK(*) or C the differential equation in subroutine DF. Any such alteration C constitutes a new problem and must be treated as such, i.e. C you must start afresh. C C You cannot change from vector to scalar error control or vice C versa (INFO(2)) but you can change the size of the entries of C RTOL, ATOL. Increasing a tolerance makes the equation easier C to integrate. Decreasing a tolerance will make the equation C harder to integrate and should generally be avoided. C C You can switch from the intermediate-output mode to the C interval mode (INFO(3)) or vice versa at any time. C C If it has been necessary to prevent the integration from going C past a point TSTOP (INFO(4), RWORK(1)), keep in mind that the C code will not integrate to any TOUT beyond the currently C specified TSTOP. Once TSTOP has been reached you must change C the value of TSTOP or set INFO(4)=0. You may change INFO(4) C or TSTOP at any time but you must supply the value of TSTOP in C RWORK(1) whenever you set INFO(4)=1. C C The parameter INFO(1) is used by the code to indicate the C beginning of a new problem and to indicate whether integration C is to be continued. you must input the value INFO(1) = 0 C when starting a new problem. You must input the value C INFO(1) = 1 if you wish to continue after an interrupted task. C Do not set INFO(1) = 0 on a continuation call unless you C want the code to restart at the current T. C C *** Following A Completed Task *** C If C IDID = 1, call the code again to continue the integration C another step in the direction of TOUT. C C IDID = 2 or 3, define a new TOUT and call the code again. C TOUT must be different from T. You cannot change C the direction of integration without restarting. C C *** Following An Interrupted Task *** C To show the code that you realize the task was C interrupted and that you want to continue, you C must take appropriate action and reset INFO(1) = 1 C If C IDID = -1, the code has attempted 500 steps. C If you want to continue, set INFO(1) = 1 and C call the code again. An additional 500 steps C will be allowed. C C IDID = -2, the error tolerances RTOL, ATOL have been C increased to values the code estimates appropriate C for continuing. You may want to change them C yourself. If you are sure you want to continue C with relaxed error tolerances, set INFO(1)=1 and C call the code again. C C IDID = -3, a solution component is zero and you set the C corresponding component of ATOL to zero. If you C are sure you want to continue, you must first C alter the error criterion to use positive values C for those components of ATOL corresponding to zero C solution components, then set INFO(1)=1 and call C the code again. C C IDID = -4, the problem appears to be stiff. It is very C inefficient to solve such problems with DDEABM. C The code DDEBDF in DEPAC handles this task C efficiently. If you are absolutely sure you want C to continue with DDEABM, set INFO(1)=1 and call C the code again. C C IDID = -5,-6,-7,..,-32 --- cannot occur with this code C but used by other members of DEPAC or possible C future extensions. C C *** Following A Terminated Task *** C If C IDID = -33, you cannot continue the solution of this C problem. An attempt to do so will result in your C run being terminated. C C ********************************************************************** C***LONG DESCRIPTION C C ********************************************************************** C * DEPAC Package Overview * C ********************************************************************** C C .... You have a choice of three differential equation solvers from C .... DEPAC. The following brief descriptions are meant to aid you in C .... choosing the most appropriate code for your problem. C C .... DDERKF is a fifth order Runge-Kutta code. It is the simplest of C .... the three choices, both algorithmically and in the use of the C .... code. DDERKF is primarily designed to solve non-stiff and C .... mildly stiff differential equations when derivative evaluations C .... are not expensive. It should generally not be used to get high C .... accuracy results nor answers at a great many specific points. C .... Because DDERKF has very low overhead costs, it will usually C .... result in the least expensive integration when solving C .... problems requiring a modest amount of accuracy and having C .... equations that are not costly to evaluate. DDERKF attempts to C .... discover when it is not suitable for the task posed. C C .... DDEABM is a variable order (one through twelve) Adams code. C .... Its complexity lies somewhere between that of DDERKF and C .... DDEBDF. DDEABM is primarily designed to solve non-stiff and C .... mildly stiff differential equations when derivative evaluations C .... are expensive, high accuracy results are needed or answers at C .... many specific points are required. DDEABM attempts to discover C .... when it is not suitable for the task posed. C C .... DDEBDF is a variable order (one through five) backward C .... differentiation formula code. it is the most complicated of C .... the three choices. DDEBDF is primarily designed to solve stiff C .... differential equations at crude to moderate tolerances. C .... If the problem is very stiff at all, DDERKF and DDEABM will be C .... quite inefficient compared to DDEBDF. However, DDEBDF will be C .... inefficient compared to DDERKF and DDEABM on non-stiff problems C .... because it uses much more storage, has a much larger overhead, C .... and the low order formulas will not give high accuracies C .... efficiently. C C .... The concept of stiffness cannot be described in a few words. C .... If you do not know the problem to be stiff, try either DDERKF C .... or DDEABM. Both of these codes will inform you of stiffness C .... when the cost of solving such problems becomes important. C C ********************************************************************* C***REFERENCES SHAMPINE, L. F., WATTS, H. A., *DEPAC - DESIGN OF A USER C ORIENTED PACKAGE OF ODE SOLVERS*, SAND79-2374, C SANDIA LABORATORIES, 1979. C***ROUTINES CALLED DDE2,XERRWV C***END PROLOGUE DDEABM C INTEGER IALPHA, IBETA, IDELSN, IDID, IFOURU, IG, IHOLD, INFLOP, 1 INFO, IP, IPAR, IPHI, IPSI, ISIG, ITOLD, ITSTAR, ITWOU, 2 IV, IW, IWORK, IWT, IYP, IYPOUT, IYY, LIW, LRW, NEQ REAL RTEMP DOUBLE PRECISION ATOL, RPAR, RTOL, RWORK, T, TOUT, Y LOGICAL START,PHASE1,NORND,STIFF,INTOUT C DIMENSION Y(NEQ),INFO(15),RTOL(1),ATOL(1),RWORK(LRW),IWORK(LIW), 1 RPAR(1),IPAR(1) C EXTERNAL DF C C .................................................................. C C INITIALIZE A COUNTER FOR MONITORING AN INFINITE LOOP PITFALL C DATA INFLOP /0/ C C .................................................................. C C CHECK FOR AN APPARENT INFINITE LOOP C C***FIRST EXECUTABLE STATEMENT DDEABM IF (INFLOP .LT. 5) GO TO 10 IF (T .NE. RWORK(21+NEQ)) GO TO 10 RTEMP = T CALL XERRWV(246HDDEABM - AN APPARENT INFINITE LOOP HAS BEEN DETEC 1TED. YOU HAVE MADE REPEATED CALLS AT T = (R1) AND INTEG 2RATION HAS NOT ADVANCED. CHECK THE WAY YOU HAVE SET PARAM 3ETERS FOR THE CALL TO THE CODE, PARTICULARLY INFO(1). 4,246,13,2,0,0,0,1,RTEMP,0.0) GO TO 50 10 CONTINUE C C CHECK LRW AND LIW FOR SUFFICIENT STORAGE ALLOCATION C IDID = 0 IF (LRW .GE. 120 + 21*NEQ) GO TO 20 CALL XERRWV(117HDDEABM - LENGTH OF RWORK ARRAY MUST BE AT LEA 1ST 120 + 21*NEQ. YOU HAVE CALLED THE CODE WITH LRW = 2 (I1).,117,1,1,1,LRW,0,0,0.0,0.0) IDID = -33 20 CONTINUE C IF (LIW .GE. 33) GO TO 30 CALL XERRWV(105HDDEABM - LENGTH OF IWORK ARRAY MUST BE AT LEAS 1T 33. YOU HAVE CALLED THE CODE WITH LIW = (I1).,105,2 2,1,1,LIW,0,0,0.0,0.0) IDID = -33 30 CONTINUE C C COMPUTE THE INDICES FOR THE ARRAYS TO BE STORED IN THE WORK C ARRAY C IYPOUT = 21 ITSTAR = NEQ + 21 IYP = 1 + ITSTAR IYY = NEQ + IYP IWT = NEQ + IYY IP = NEQ + IWT IPHI = NEQ + IP IALPHA = (NEQ*16) + IPHI IBETA = 12 + IALPHA IPSI = 12 + IBETA IV = 12 + IPSI IW = 12 + IV ISIG = 12 + IW IG = 13 + ISIG IHOLD = 13 + IG ITOLD = 1 + IHOLD IDELSN = 1 + ITOLD ITWOU = 1 + IDELSN IFOURU = 1 + ITWOU C RWORK(ITSTAR) = T IF (INFO(1) .EQ. 0) GO TO 40 START = IWORK(21) .NE. (-1) PHASE1 = IWORK(22) .NE. (-1) NORND = IWORK(23) .NE. (-1) STIFF = IWORK(24) .NE. (-1) INTOUT = IWORK(25) .NE. (-1) 40 CONTINUE C CALL DDE2(DF,NEQ,T,Y,TOUT,INFO,RTOL,ATOL,IDID,RWORK(IYPOUT), 1 RWORK(IYP),RWORK(IYY),RWORK(IWT),RWORK(IP), 2 RWORK(IPHI),RWORK(IALPHA),RWORK(IBETA),RWORK(IPSI), 3 RWORK(IV),RWORK(IW),RWORK(ISIG),RWORK(IG),RWORK(11), 4 RWORK(12),RWORK(13),RWORK(IHOLD),RWORK(ITOLD), 5 RWORK(IDELSN),RWORK(1),RWORK(ITWOU),RWORK(IFOURU), 6 START,PHASE1,NORND,STIFF,INTOUT,IWORK(26),IWORK(27), 7 IWORK(28),IWORK(29),IWORK(30),IWORK(31),IWORK(32), 8 RPAR,IPAR) C IWORK(21) = -1 IF (START) IWORK(21) = 1 IWORK(22) = -1 IF (PHASE1) IWORK(22) = 1 IWORK(23) = -1 IF (NORND) IWORK(23) = 1 IWORK(24) = -1 IF (STIFF) IWORK(24) = 1 IWORK(25) = -1 IF (INTOUT) IWORK(25) = 1 C IF (IDID .NE. (-2)) INFLOP = INFLOP + 1 IF (T .NE. RWORK(ITSTAR)) INFLOP = 0 50 CONTINUE C RETURN END SUBROUTINE DDE2(DF,NEQ,T,Y,TOUT,INFO,RTOL,ATOL,IDID,YPOUT,YP,YY, 1 WT,P,PHI,ALPHA,BETA,PSI,V,W,SIG,G,H,EPS,X,HOLD,TOLD,DELSGN, 2 TSTOP,TWOU,FOURU,START,PHASE1,NORND,STIFF,INTOUT,NS,KORD,KOLD, 3 INIT,KSTEPS,KLE4,IQUIT,RPAR,IPAR) C***BEGIN PROLOGUE DDE2 C***REFER TO DDEABM C C DDEABM merely allocates storage for DDE2 to relieve the user of C the inconvenience of a long call list. Consequently DDE2 is used C as described in the comments for DDEABM . C***ROUTINES CALLED D1MACH,DINTRP,DSTEP2,XERRWV C***END PROLOGUE DDE2 C INTEGER IDID, INFO, INIT, IPAR, IQUIT, K, KLE4, KOLD, KORD, 1 KSTEPS, L, LTOL, MAXNUM, NATOLP, NEQ, NRTOLP, NS REAL RTEMP1, RTEMP2 DOUBLE PRECISION A, ABSDEL, ALPHA, ATOL, BETA, D1MACH, DABS, 1 DEL, DELSGN, DMAX1, DMIN1, DSIGN, DT, EPS, FOURU, G, H, 2 HA, HOLD, P, PHI, PSI, RPAR, RTOL, SIG, T, TOLD, TOUT, 3 TSTOP, TWOU, U, V, W, WT, X, Y, YP, YPOUT, YY LOGICAL STIFF,CRASH,START,PHASE1,NORND,INTOUT C DIMENSION Y(NEQ),YY(NEQ),WT(NEQ),PHI(NEQ,16),P(NEQ),YP(NEQ), 1 YPOUT(NEQ),PSI(12),ALPHA(12),BETA(12),SIG(13),V(12), 2 W(12),G(13),INFO(15),RTOL(1),ATOL(1),RPAR(1),IPAR(1) C EXTERNAL DF C C .................................................................. C C THE EXPENSE OF SOLVING THE PROBLEM IS MONITORED BY COUNTING THE C NUMBER OF STEPS ATTEMPTED. WHEN THIS EXCEEDS MAXNUM, THE C COUNTER IS RESET TO ZERO AND THE USER IS INFORMED ABOUT POSSIBLE C EXCESSIVE WORK. C DATA MAXNUM /500/ C C .................................................................. C C***FIRST EXECUTABLE STATEMENT DDE2 IF (INFO(1) .NE. 0) GO TO 10 C C ON THE FIRST CALL , PERFORM INITIALIZATION -- C DEFINE THE MACHINE UNIT ROUNDOFF QUANTITY U BY CALLING C THE FUNCTION ROUTINE D1MACH. THE USER MUST MAKE SURE C THAT THE VALUES SET IN D1MACH ARE RELEVANT TO THE C COMPUTER BEING USED. C U = D1MACH(4) C -- SET ASSOCIATED MACHINE DEPENDENT PARAMETERS TWOU = 2.0D0*U FOURU = 4.0D0*U C -- SET TERMINATION FLAG IQUIT = 0 C -- SET INITIALIZATION INDICATOR INIT = 0 C -- SET COUNTER FOR ATTEMPTED STEPS KSTEPS = 0 C -- SET INDICATOR FOR INTERMEDIATE-OUTPUT INTOUT = .FALSE. C -- SET INDICATOR FOR STIFFNESS DETECTION STIFF = .FALSE. C -- SET STEP COUNTER FOR STIFFNESS DETECTION KLE4 = 0 C -- SET INDICATORS FOR DSTEP2 CODE START = .TRUE. PHASE1 = .TRUE. NORND = .TRUE. C -- RESET INFO(1) FOR SUBSEQUENT CALLS INFO(1) = 1 10 CONTINUE C C .................................................................. C C CHECK VALIDITY OF INPUT PARAMATERS ON EACH ENTRY C IF (INFO(1) .EQ. 0 .OR. INFO(1) .EQ. 1) GO TO 20 CALL XERRWV(249HDDEABM - INFO(1) MUST BE SET TO 0 FOR THE START O 1F A NEW PROBLEM, AND MUST BE SET TO 1 FOLLOWING AN INTERR 2UPTED TASK. YOU ARE ATTEMPTING TO CONTINUE THE INTEG 3RATION ILLEGALLY BY CALLING THE CODE WITH INFO(1) = (I1 4).,249,3, 1,1,INFO(1),0,0,0.0,0.0) IDID = -33 20 CONTINUE C IF (INFO(2) .EQ. 0 .OR. INFO(2) .EQ. 1) GO TO 30 CALL XERRWV(164HDDEABM - INFO(2) MUST BE 0 OR 1 INDICATING SCALA 1R AND VECTOR ERROR TOLERANCES, RESPECTIVELY. YOU HAVE 2 CALLED THE CODE WITH INFO(2) = (I1).,164,4,1,1,INFO 3(2),0,0, 0.0,0.0) IDID = -33 30 CONTINUE C IF (INFO(3) .EQ. 0 .OR. INFO(3) .EQ. 1) GO TO 40 CALL XERRWV(195HDDEABM - INFO(3) MUST BE 0 OR 1 INDICATING THE I 1NTERVAL OR INTERMEDIATE-OUTPUT MODE OF INTEGRAT 2ION, RESPECTIVELY. YOU HAVE CALLED THE CODE WITH I 3NFO(3) = (I1).,195,5,1,1,INFO(3),0,0,0.0,0.0) IDID = -33 40 CONTINUE C IF (INFO(4) .EQ. 0 .OR. INFO(4) .EQ. 1) GO TO 50 CALL XERRWV(195HDDEABM - INFO(4) MUST BE 0 OR 1 INDICATING WHETHE 1R OR NOT THE INTEGRATION INTERVAL IS TO BE RESTR 2ICTED BY A POINT TSTOP. YOU HAVE CALLED THE CODE WITH I 3NFO(4) = (I1).,195,14,1,1,INFO(4),0,0,0.0,0.0) IDID = -33 50 CONTINUE C IF (NEQ .GE. 1) GO TO 60 CALL XERRWV(117HDDEABM - THE NUMBER OF EQUATIONS NEQ MUST BE A PO 1SITIVE INTEGER. YOU HAVE CALLED THE CODE WITH NEQ = (I1 2).,117,6, 1,1,NEQ,0,0,0.0,0.0) IDID = -33 60 CONTINUE C NRTOLP = 0 NATOLP = 0 DO 110 K = 1, NEQ C BEGIN BLOCK PERMITTING ...EXITS TO 100 IF (NRTOLP .GT. 0) GO TO 90 IF (RTOL(K) .GE. 0.0D0) GO TO 70 RTEMP1 = RTOL(K) CALL XERRWV(238HDDEABM - THE RELATIVE ERROR TOLERANCES 1 RTOL MUST BE NON-NEGATIVE. YOU HAVE CALLED THE CODE WITH 2 RTOL(I1) = (R1). IN THE CASE OF VECTOR ERROR TOLERANC 3ES, NO FURTHER CHECKING OF RTOL COMPONENTS IS DON 4E.,238,7,1, 1,K,0,1,RTEMP1,0.0) IDID = -33 C ............EXIT IF (NATOLP .GT. 0) GO TO 120 NRTOLP = 1 GO TO 80 70 CONTINUE C .........EXIT IF (NATOLP .GT. 0) GO TO 100 80 CONTINUE 90 CONTINUE C ...EXIT IF (ATOL(K) .GE. 0.0D0) GO TO 100 RTEMP1 = ATOL(K) CALL XERRWV(238HDDEABM - THE ABSOLUTE ERROR TOLERANCES ATOL MU 1ST BE NON-NEGATIVE. YOU HAVE CALLED THE CODE WITH ATOL(I 21) = (R1). IN THE CASE OF VECTOR ERROR TOLERANCES, NO F 3URTHER CHECKING OF ATOL COMPONENTS IS DONE.,238, 48,1,1,K,0, 1,RTEMP1,0.0) IDID = -33 C ......EXIT IF (NRTOLP .GT. 0) GO TO 120 NATOLP = 1 100 CONTINUE C ...EXIT IF (INFO(2) .EQ. 0) GO TO 120 110 CONTINUE 120 CONTINUE C IF (INFO(4) .NE. 1) GO TO 130 IF (DSIGN(1.0D0,TOUT-T) .EQ. DSIGN(1.0D0,TSTOP-T) 1 .AND. DABS(TOUT-T) .LE. DABS(TSTOP-T)) GO TO 130 RTEMP1 = TOUT RTEMP2 = TSTOP CALL XERRWV(191HDDEABM - YOU HAVE CALLED THE CODE WITH TOUT = (R 11) BUT YOU HAVE ALSO TOLD THE CODE (INFO(4) = 1) NOT TO I 2NTEGRATE PAST THE POINT TSTOP = (R2). THESE INSTRUCTIONS 3 CONFLICT.,191,14,1,0,0,0,2,RTEMP1,RTEMP2) IDID = -33 130 CONTINUE C IF (INIT .EQ. 0) GO TO 170 C CHECK SOME CONTINUATION POSSIBILITIES IF (T .NE. TOUT) GO TO 140 RTEMP1 = T CALL XERRWV(116HDDEABM - YOU HAVE CALLED THE CODE WITH T = TO 1UT AT T = (R1). THIS IS NOT ALLOWED ON CONTINUATION CA 2LLS.,116,9,1,0,0,0,1,RTEMP1,0.0) IDID = -33 140 CONTINUE C IF (T .EQ. TOLD) GO TO 150 RTEMP1 = TOLD RTEMP2 = T CALL XERRWV(113HDDEABM - YOU HAVE CHANGED THE VALUE OF T FROM 1 (R1) TO (R2). THIS IS NOT ALLOWED ON CONTINUATION CALL 2S.,113,10, 1,0,0,0,2,RTEMP1,RTEMP2) IDID = -33 150 CONTINUE C IF (INIT .EQ. 1) GO TO 160 IF (DELSGN*(TOUT - T) .GE. 0.0D0) GO TO 160 RTEMP1 = TOUT CALL XERRWV(168HDDEABM - BY CALLING THE CODE WITH TOUT = (R1) 1 , YOU ARE ATTEMPTING TO CHANGE THE DIRECTION OF INTEGRAT 2ION. THIS IS NOT ALLOWED WITHOUT RESTARTING.,168,11,1 3,0,0,0,1, RTEMP1,0.0) IDID = -33 160 CONTINUE 170 CONTINUE C IF (IDID .NE. (-33)) GO TO 200 IF (IQUIT .EQ. (-33)) GO TO 180 C C INVALID INPUT DETECTED IQUIT = -33 INFO(1) = -1 GO TO 190 180 CONTINUE C CALL XERRWV(191HDDEABM - INVALID INPUT WAS DETECTED ON SUCCESS 1IVE ENTRIES. IT IS IMPOSSIBLE TO PROCEED BECAUSE YO 2U HAVE NOT CORRECTED THE PROBLEM, SO EXECUTION IS BEIN 3G TERMINATED.,191,12,2,0,0,0,0,0.0,0.0) 190 CONTINUE GO TO 500 200 CONTINUE C C ............................................................... C C RTOL = ATOL = 0. IS ALLOWED AS VALID INPUT AND INTERPRETED C AS ASKING FOR THE MOST ACCURATE SOLUTION POSSIBLE. IN THIS C CASE, THE RELATIVE ERROR TOLERANCE RTOL IS RESET TO THE C SMALLEST VALUE FOURU WHICH IS LIKELY TO BE REASONABLE FOR C THIS METHOD AND MACHINE C DO 220 K = 1, NEQ IF (RTOL(K) + ATOL(K) .GT. 0.0D0) GO TO 210 RTOL(K) = FOURU IDID = -2 210 CONTINUE C ...EXIT IF (INFO(2) .EQ. 0) GO TO 230 220 CONTINUE 230 CONTINUE C IF (IDID .NE. (-2)) GO TO 240 C RTOL=ATOL=0 ON INPUT, SO RTOL IS CHANGED TO A C SMALL POSITIVE VALUE INFO(1) = -1 GO TO 490 240 CONTINUE C BEGIN BLOCK PERMITTING ...EXITS TO 480 C BEGIN BLOCK PERMITTING ...EXITS TO 290 C BEGIN BLOCK PERMITTING ...EXITS TO 270 C C BRANCH ON STATUS OF INITIALIZATION INDICATOR C INIT=0 MEANS INITIAL DERIVATIVES AND NOMINAL C STEP SIZE C AND DIRECTION NOT YET SET C INIT=1 MEANS NOMINAL STEP SIZE AND DIRECTION C NOT YET SET INIT=2 MEANS NO FUTHER C INITIALIZATION REQUIRED C IF (INIT .EQ. 0) GO TO 250 C ......EXIT IF (INIT .EQ. 1) GO TO 270 C .........EXIT GO TO 290 250 CONTINUE C C ................................................... C C MORE INITIALIZATION -- C -- EVALUATE INITIAL C DERIVATIVES C INIT = 1 A = T CALL DF(A,Y,YP,RPAR,IPAR) C ...EXIT IF (T .NE. TOUT) GO TO 270 IDID = 2 DO 260 L = 1, NEQ YPOUT(L) = YP(L) 260 CONTINUE TOLD = T C .........EXIT GO TO 480 270 CONTINUE C C -- SET INDEPENDENT AND DEPENDENT VARIABLES C X AND YY(*) FOR DSTEP2 C -- SET DSIGN OF INTEGRATION DIRECTION C -- INITIALIZE THE STEP SIZE C INIT = 2 X = T DO 280 L = 1, NEQ YY(L) = Y(L) 280 CONTINUE DELSGN = DSIGN(1.0D0,TOUT-T) H = DSIGN(DMAX1(FOURU*DABS(X),DABS(TOUT-X)),TOUT-X) 290 CONTINUE C C ......................................................... C C ON EACH CALL SET INFORMATION WHICH DETERMINES THE C ALLOWED INTERVAL OF INTEGRATION BEFORE RETURNING WITH C AN ANSWER AT TOUT C DEL = TOUT - T ABSDEL = DABS(DEL) C C ......................................................... C C IF ALREADY PAST OUTPUT POINT, INTERPOLATE AND RETURN C 300 CONTINUE C BEGIN BLOCK PERMITTING ...EXITS TO 430 C BEGIN BLOCK PERMITTING ...EXITS TO 410 IF (DABS(X-T) .LT. ABSDEL) GO TO 320 CALL DINTRP(X,YY,TOUT,Y,YPOUT,NEQ,KOLD,PHI, 1 PSI) IDID = 3 IF (X .NE. TOUT) GO TO 310 IDID = 2 INTOUT = .FALSE. 310 CONTINUE T = TOUT TOLD = T C ...............EXIT GO TO 480 320 CONTINUE C C IF CANNOT GO PAST TSTOP AND SUFFICIENTLY CLOSE, C EXTRAPOLATE AND RETURN C IF (INFO(4) .NE. 1) GO TO 340 IF (DABS(TSTOP-X) .GE. FOURU*DABS(X)) GO TO 340 DT = TOUT - X DO 330 L = 1, NEQ Y(L) = YY(L) + DT*YP(L) 330 CONTINUE CALL DF(TOUT,Y,YPOUT,RPAR,IPAR) IDID = 3 T = TOUT TOLD = T C ...............EXIT GO TO 480 340 CONTINUE C IF (INFO(3) .EQ. 0 .OR. .NOT.INTOUT) GO TO 360 C C INTERMEDIATE-OUTPUT MODE C IDID = 1 DO 350 L = 1, NEQ Y(L) = YY(L) YPOUT(L) = YP(L) 350 CONTINUE T = X TOLD = T INTOUT = .FALSE. C ...............EXIT GO TO 480 360 CONTINUE C C ................................................ C C MONITOR NUMBER OF STEPS ATTEMPTED C IF (KSTEPS .LE. MAXNUM) GO TO 390 C C A SIGNIFICANT AMOUNT OF WORK HAS BEEN C EXPENDED IDID = -1 KSTEPS = 0 IF (.NOT.STIFF) GO TO 370 C C PROBLEM APPEARS TO BE STIFF IDID = -4 STIFF = .FALSE. KLE4 = 0 370 CONTINUE C DO 380 L = 1, NEQ Y(L) = YY(L) YPOUT(L) = YP(L) 380 CONTINUE T = X TOLD = T INFO(1) = -1 INTOUT = .FALSE. C ...............EXIT GO TO 480 390 CONTINUE C C ................................................ C C LIMIT STEP SIZE, SET WEIGHT VECTOR AND TAKE A C STEP C HA = DABS(H) IF (INFO(4) .EQ. 1) HA = DMIN1(HA,DABS(TSTOP-X)) H = DSIGN(HA,H) EPS = 1.0D0 LTOL = 1 DO 400 L = 1, NEQ IF (INFO(2) .EQ. 1) LTOL = L WT(L) = RTOL(LTOL)*DABS(YY(L)) + ATOL(LTOL) C ......EXIT IF (WT(L) .LE. 0.0D0) GO TO 410 400 CONTINUE C ......EXIT GO TO 430 410 CONTINUE C C RELATIVE ERROR CRITERION INAPPROPRIATE IDID = -3 DO 420 L = 1, NEQ Y(L) = YY(L) YPOUT(L) = YP(L) 420 CONTINUE T = X TOLD = T INFO(1) = -1 INTOUT = .FALSE. C .........EXIT GO TO 480 430 CONTINUE C CALL DSTEP2(DF,NEQ,YY,X,H,EPS,WT,START,HOLD,KORD,KOLD, 1 CRASH,PHI,P,YP,PSI,ALPHA,BETA,SIG,V,W,G, 2 PHASE1,NS,NORND,KSTEPS,TWOU,FOURU,RPAR, 3 IPAR) C C ...................................................... C IF (.NOT.CRASH) GO TO 470 C C TOLERANCES TOO SMALL IDID = -2 RTOL(1) = EPS*RTOL(1) ATOL(1) = EPS*ATOL(1) IF (INFO(2) .EQ. 0) GO TO 450 DO 440 L = 2, NEQ RTOL(L) = EPS*RTOL(L) ATOL(L) = EPS*ATOL(L) 440 CONTINUE 450 CONTINUE DO 460 L = 1, NEQ Y(L) = YY(L) YPOUT(L) = YP(L) 460 CONTINUE T = X TOLD = T INFO(1) = -1 INTOUT = .FALSE. C .........EXIT GO TO 480 470 CONTINUE C C (STIFFNESS TEST) COUNT NUMBER OF CONSECUTIVE STEPS C TAKEN WITH THE ORDER OF THE METHOD BEING LESS OR EQUAL C TO FOUR C KLE4 = KLE4 + 1 IF (KOLD .GT. 4) KLE4 = 0 IF (KLE4 .GE. 50) STIFF = .TRUE. INTOUT = .TRUE. GO TO 300 480 CONTINUE 490 CONTINUE 500 CONTINUE RETURN END SUBROUTINE DSTEP2(DF,NEQN,Y,X,H,EPS,WT,START,HOLD,K,KOLD,CRASH, 1 PHI,P,YP,PSI,ALPHA,BETA,SIG,V,W,G,PHASE1,NS,NORND,KSTEPS,TWOU, 2 FOURU,RPAR,IPAR) C***BEGIN PROLOGUE DSTEP2 C***DATE WRITTEN 820301 (YYMMDD) C***REVISION DATE 830701 (YYMMDD) C***CATEGORY NO. I1A1B C***KEYWORDS ADAMS METHOD,DEPAC,INITIAL VALUE PROBLEMS,ODE, C ORDINARY DIFFERENTIAL EQUATIONS,PREDICTOR-CORRECTOR C***AUTHOR SHAMPINE, L. F., (SNLA) C GORDON, M. K. C***PURPOSE Integrates a system of first order ODES one step. C***DESCRIPTION C C Written by L. F. Shampine and M. K. Gordon C C Abstract C C Subroutine DSTEP2 is normally used indirectly through subroutine C DDEABM . Because DDEABM suffices for most problems and is much C easier to use, using it should be considered before using DSTEP2 C alone. C C Subroutine DSTEP2 integrates a system of NEQN first order ordinary C differential equations one step, normally from X to X+H, using a C modified divided difference form of the Adams Pece formulas. Local C extrapolation is used to improve absolute stability and accuracy. C The code adjusts its order and step size to control the local error C per unit step in a generalized sense. Special devices are included C to control roundoff error and to detect when the user is requesting C too much accuracy. C C This code is completely explained and documented in the text, C COMPUTER SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS, THE INITIAL C VALUE PROBLEM by L. F. Shampine and M. K. Gordon. C Further details on use of this code are available in *SOLVING C ORDINARY DIFFERENTIAL EQUATIONS WITH ODE, STEP, AND DINTRP*, C by L. F. Shampine and M. K. Gordon, SLA-73-1060. C C C The parameters represent -- C DF -- Subroutine to evaluate derivatives C NEQN -- Number of equations to be integrated C Y(*) -- Solution vector at X C X -- Independent variable C H -- Appropriate step size for next step. Normally determined by C code C EPS -- Local error tolerance C WT(*) -- Vector of weights for error criterion C START -- LOGICAL variable set .TRUE. for first step, .FALSE. C otherwise C HOLD -- Step size used for last successful step C K -- Appropriate order for next step (determined by code) C KOLD -- Order used for last successful step C CRASH -- LOGICAL variable set .TRUE. when no step can be taken, C .FALSE. otherwise. C YP(*) -- Derivative of solution vector at X after successful C step C KSTEPS -- Counter on attempted steps C TWOU -- 2.*U where U is machine unit roundoff quantity C FOURU -- 4.*U where U is machine unit roundoff quantity C RPAR,IPAR -- Parameter arrays which you may choose to use C for communication between your program and subroutine DF. C They are not altered or used by DSTEP2. C The arrays PHI, PSI are required for the interpolation subroutine C DINTRP . The array P is internal to the code. The remaining nine C variables and arrays are included in the call list only to eliminate C local retention of variables between calls. C C Input to DSTEP2 C C First call -- C C The user must provide storage in his calling program for all arrays C in the call list, namely C C DIMENSION Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),PSI(12), C 1 ALPHA(12),BETA(12),SIG(13),V(12),W(12),G(13),RPAR(1),IPAR(1) C -- -- **NOTE** C C The user must also declare START , CRASH , PHASE1 and NORND C LOGICAL variables and DF an external subroutine, supply the C Subroutine DF(X,Y,YP) to evaluate C DY(I)/DX = YP(I) = DF(X,Y(1),Y(2),...,Y(NEQN)) C and initialize only the following parameters. C NEQN -- Number of equations to be integrated C Y(*) -- Vector of initial values of dependent variables C X -- Initial value of the independent variable C H -- Nominal step size indicating direction of integration C and maximum size of step. Must be variable C EPS -- Local error tolerance per step. Must be variable C WT(*) -- Vector of non-zero weights for error criterion C START -- .TRUE. C KSTEPS -- Set KSTEPS to zero C TWOU -- 2.*U where U is machine unit roundoff quantity C FOURU -- 4.*U where U is machine unit roundoff quantity C Define U to be the machine unit roundoff quantity by calling C the function routine D1MACH, U = D1MACH(4), or by C computing U so that u is the smallest positive number such C that 1.0+U .GT. 1.0. C C DSTEP2 requires that the L2 norm of the vector with components C local ERROR(L)/WT(L) be less than EPS for a successful step. The C array WT allows the user to specify an error test appropriate C for his problem. For example, C WT(L) = 1.0 specifies absolute error, C = DABS(Y(L)) error relative to the most recent value of C the C L-th component of the solution, C = DABS(YP(L)) error relative to the most recent value of C the L-th component of the derivative, C = DMAX1(WT(L),DABS(Y(L))) error relative to the largest C magnitude of L-th component obtained so far, C = DABS(Y(L))*RELERR/EPS + ABSERR/EPS specifies a mixed C relative-absolute test where RELERR is relative C error, ABSERR is absolute error and EPS = C DMAX1(RELERR,ABSERR) . C C Subsequent calls -- C C Subroutine DSTEP2 is designed so that all information needed to C continue the integration, including the step size H and the order C K , is returned with each step. With the exception of the step C size, the error tolerance, and the weights, none of the parameters C should be altered. The array WT must be updated after each step C to maintain relative error tests like those above. Normally the C integration is continued just beyond the desired endpoint and the C solution interpolated there with subroutine DINTRP . If it is C impossible to integrate beyond the endpoint, the step size may be C reduced to hit the endpoint since the code will not take a step C larger than the H input. Changing the direction of integration, C i.e., the DSIGN of H , requires the user set START = .TRUE. before C calling DSTEP2 again. This is the only situation in which START C should be altered. C C Output from DSTEP2 C C Successful step -- C C The subroutine returns after each successful step with START and C CRASH set .FALSE. . X represents the independent variable C advanced one step of length HOLD from its value on input and Y C the solution vector at the new value of X . All other parameters C represent information corresponding to the new X needed to C continue the integration. C C Unsuccessful step -- C C When the error tolerance is too small for the machine precision, C the subroutine returns without taking a step and CRASH = .TRUE. . C An appropriate step size and error tolerance for continuing are C estimated and all other information is restored as upon input C before returning. To continue with the larger tolerance, the user C just calls the code again. A restart is neither required nor C desirable. C***REFERENCES SHAMPINE, L. F., GORDON, M. K., *SOLVING ORDINARY C DIFFERENTIAL EQUATIONS WITH ODE, STEP, AND DINTRP*, C SLA-73-1060, SANDIA LABORATORIES, 1973. C***ROUTINES CALLED D1MACH,DHSTRT C***END PROLOGUE DSTEP2 C INTEGER I, IFAIL, IM1, IP1, IPAR, IQ, J, K, KM1, KM2, KNEW, 1 KOLD, KP1, KP2, KSTEPS, L, LIMIT1, LIMIT2, NEQN, NS, NSM2, 2 NSP1, NSP2 DOUBLE PRECISION ABSH, ALPHA, BETA, BIG, D1MACH, DABS, DMAX1, 1 DMIN1, DSIGN, DSQRT, EPS, ERK, ERKM1, ERKM2, ERKP1, ERR, 2 FOURU, G, GSTR, H, HNEW, HOLD, P, P5EPS, PHI, PSI, R, 3 REALI, REALNS, RHO, ROUND, RPAR, SIG, TAU, TEMP1, TEMP2, 4 TEMP3, TEMP4, TEMP5, TEMP6, TWO, TWOU, U, V, W, WT, X, 5 XOLD, Y, YP LOGICAL START,CRASH,PHASE1,NORND DIMENSION Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),PSI(12), 1 ALPHA(12),BETA(12),SIG(13),V(12),W(12),G(13),RPAR(1), 2 IPAR(1) DIMENSION TWO(13),GSTR(13) EXTERNAL DF C DATA TWO(1),TWO(2),TWO(3),TWO(4),TWO(5),TWO(6),TWO(7),TWO(8), 1 TWO(9),TWO(10),TWO(11),TWO(12),TWO(13) 2 /2.0D0,4.0D0,8.0D0,16.0D0,32.0D0,64.0D0,128.0D0,256.0D0, 3 512.0D0,1024.0D0,2048.0D0,4096.0D0,8192.0D0/ DATA GSTR(1),GSTR(2),GSTR(3),GSTR(4),GSTR(5),GSTR(6),GSTR(7), 1 GSTR(8),GSTR(9),GSTR(10),GSTR(11),GSTR(12),GSTR(13) 2 /0.5D0,0.0833D0,0.0417D0,0.0264D0,0.0188D0,0.0143D0,0.0114D0, 3 0.00936D0,0.00789D0,0.00679D0,0.00592D0,0.00524D0,0.00468D0/ C C C *** BEGIN BLOCK 0 *** C CHECK IF STEP SIZE OR ERROR TOLERANCE IS TOO SMALL FOR MACHINE C PRECISION. IF FIRST STEP, INITIALIZE PHI ARRAY AND ESTIMATE A C STARTING STEP SIZE. C *** C C IF STEP SIZE IS TOO SMALL, DETERMINE AN ACCEPTABLE ONE C C***FIRST EXECUTABLE STATEMENT DSTEP2 CRASH = .TRUE. IF (DABS(H) .GE. FOURU*DABS(X)) GO TO 10 H = DSIGN(FOURU*DABS(X),H) GO TO 650 10 CONTINUE P5EPS = 0.5D0*EPS C C IF ERROR TOLERANCE IS TOO SMALL, INCREASE IT TO AN ACCEPTABLE C VALUE C ROUND = 0.0D0 DO 20 L = 1, NEQN ROUND = ROUND + (Y(L)/WT(L))**2 20 CONTINUE ROUND = TWOU*DSQRT(ROUND) IF (P5EPS .GE. ROUND) GO TO 30 EPS = 2.0D0*ROUND*(1.0D0 + FOURU) GO TO 640 30 CONTINUE C BEGIN BLOCK PERMITTING ...EXITS TO 630 C BEGIN BLOCK PERMITTING ...EXITS TO 600 C BEGIN BLOCK PERMITTING ...EXITS TO 60 CRASH = .FALSE. G(1) = 1.0D0 G(2) = 0.5D0 SIG(1) = 1.0D0 C ...EXIT IF (.NOT.START) GO TO 60 C C INITIALIZE. COMPUTE APPROPRIATE STEP SIZE FOR C FIRST STEP C CALL DF(X,Y,YP,RPAR,IPAR) C SUM = 0.0 DO 40 L = 1, NEQN PHI(L,1) = YP(L) PHI(L,2) = 0.0D0 40 CONTINUE C 20 SUM = SUM + (YP(L)/WT(L))**2 C SUM = DSQRT(SUM) C ABSH = DABS(H) C IF(EPS .LT. 16.0*SUM*H*H) ABSH = C 0.25*DSQRT(EPS/SUM) H = C DSIGN(DMAX1(ABSH,FOURU*DABS(X)),H) C U = D1MACH(4) BIG = DSQRT(D1MACH(2)) CALL DHSTRT(DF,NEQN,X,X+H,Y,YP,WT,1,U,BIG,PHI(1,3), 1 PHI(1,4),PHI(1,5),PHI(1,6),RPAR,IPAR,H) C HOLD = 0.0D0 K = 1 KOLD = 0 START = .FALSE. PHASE1 = .TRUE. NORND = .TRUE. C ...EXIT IF (P5EPS .GT. 100.0D0*ROUND) GO TO 60 NORND = .FALSE. DO 50 L = 1, NEQN PHI(L,15) = 0.0D0 50 CONTINUE 60 CONTINUE IFAIL = 0 C *** END BLOCK 0 *** C C *** BEGIN BLOCK 1 *** C COMPUTE COEFFICIENTS OF FORMULAS FOR THIS STEP. AVOID C COMPUTING THOSE QUANTITIES NOT CHANGED WHEN STEP SIZE C IS NOT CHANGED. C *** C 70 CONTINUE C BEGIN BLOCK PERMITTING ...EXITS TO 580 C BEGIN BLOCK PERMITTING ...EXITS TO 350 KP1 = K + 1 KP2 = K + 2 KM1 = K - 1 KM2 = K - 2 C C NS IS THE NUMBER OF STEPS TAKEN WITH SIZE H, C INCLUDING THE CURRENT ONE. WHEN K.LT.NS, NO C COEFFICIENTS CHANGE C IF (H .NE. HOLD) NS = 0 IF (NS .LE. KOLD) NS = NS + 1 NSP1 = NS + 1 IF (K .LT. NS) GO TO 200 C C COMPUTE THOSE COMPONENTS OF C ALPHA(*),BETA(*),PSI(*),SIG(*) WHICH ARE C CHANGED C BETA(NS) = 1.0D0 REALNS = NS ALPHA(NS) = 1.0D0/REALNS TEMP1 = H*REALNS SIG(NSP1) = 1.0D0 IF (K .LT. NSP1) GO TO 90 DO 80 I = NSP1, K IM1 = I - 1 TEMP2 = PSI(IM1) PSI(IM1) = TEMP1 BETA(I) = BETA(IM1)*PSI(IM1)/TEMP2 TEMP1 = TEMP2 + H ALPHA(I) = H/TEMP1 REALI = I SIG(I+1) = REALI*ALPHA(I)*SIG(I) 80 CONTINUE 90 CONTINUE PSI(K) = TEMP1 C C COMPUTE COEFFICIENTS G(*) C C INITIALIZE V(*) AND SET W(*). C IF (NS .GT. 1) GO TO 110 DO 100 IQ = 1, K TEMP3 = IQ*(IQ + 1) V(IQ) = 1.0D0/TEMP3 W(IQ) = V(IQ) 100 CONTINUE GO TO 160 110 CONTINUE C C IF ORDER WAS RAISED, UPDATE DIAGONAL C PART OF V(*) C IF (K .LE. KOLD) GO TO 140 TEMP4 = K*KP1 V(K) = 1.0D0/TEMP4 NSM2 = NS - 2 IF (NSM2 .LT. 1) GO TO 130 DO 120 J = 1, NSM2 I = K - J V(I) = V(I) - ALPHA(J+1)*V(I+1) 120 CONTINUE 130 CONTINUE 140 CONTINUE C C UPDATE V(*) AND SET W(*) C LIMIT1 = KP1 - NS TEMP5 = ALPHA(NS) DO 150 IQ = 1, LIMIT1 V(IQ) = V(IQ) - TEMP5*V(IQ+1) W(IQ) = V(IQ) 150 CONTINUE G(NSP1) = W(1) 160 CONTINUE C C COMPUTE THE G(*) IN THE WORK VECTOR W(*) C NSP2 = NS + 2 IF (KP1 .LT. NSP2) GO TO 190 DO 180 I = NSP2, KP1 LIMIT2 = KP2 - I TEMP6 = ALPHA(I-1) DO 170 IQ = 1, LIMIT2 W(IQ) = W(IQ) - TEMP6*W(IQ+1) 170 CONTINUE G(I) = W(1) 180 CONTINUE 190 CONTINUE 200 CONTINUE C *** END BLOCK 1 *** C C *** BEGIN BLOCK 2 *** C PREDICT A SOLUTION P(*), EVALUATE DERIVATIVES C USING PREDICTED SOLUTION, ESTIMATE LOCAL C ERROR AT ORDER K AND ERRORS AT ORDERS K, C K-1, K-2 AS IF CONSTANT STEP SIZE WERE USED. C *** C C INCREMENT COUNTER ON ATTEMPTED STEPS C KSTEPS = KSTEPS + 1 C C CHANGE PHI TO PHI STAR C IF (K .LT. NSP1) GO TO 230 DO 220 I = NSP1, K TEMP1 = BETA(I) DO 210 L = 1, NEQN PHI(L,I) = TEMP1*PHI(L,I) 210 CONTINUE 220 CONTINUE 230 CONTINUE C C PREDICT SOLUTION AND DIFFERENCES C DO 240 L = 1, NEQN PHI(L,KP2) = PHI(L,KP1) PHI(L,KP1) = 0.0D0 P(L) = 0.0D0 240 CONTINUE DO 260 J = 1, K I = KP1 - J IP1 = I + 1 TEMP2 = G(I) DO 250 L = 1, NEQN P(L) = P(L) + TEMP2*PHI(L,I) PHI(L,I) = PHI(L,I) + PHI(L,IP1) 250 CONTINUE 260 CONTINUE IF (NORND) GO TO 280 DO 270 L = 1, NEQN TAU = H*P(L) - PHI(L,15) P(L) = Y(L) + TAU PHI(L,16) = (P(L) - Y(L)) - TAU 270 CONTINUE GO TO 300 280 CONTINUE DO 290 L = 1, NEQN P(L) = Y(L) + H*P(L) 290 CONTINUE 300 CONTINUE XOLD = X X = X + H ABSH = DABS(H) CALL DF(X,P,YP,RPAR,IPAR) C C ESTIMATE ERRORS AT ORDERS K,K-1,K-2 C ERKM2 = 0.0D0 ERKM1 = 0.0D0 ERK = 0.0D0 DO 330 L = 1, NEQN C BEGIN BLOCK PERMITTING ...EXITS TO 320 TEMP3 = 1.0D0/WT(L) TEMP4 = YP(L) - PHI(L,1) IF (KM2 .EQ. 0) GO TO 310 C ......EXIT IF (KM2 .LT. 0) GO TO 320 ERKM2 = ERKM2 1 + ((PHI(L,KM1) + TEMP4) 2 *TEMP3)**2 310 CONTINUE ERKM1 = ERKM1 1 + ((PHI(L,K) + TEMP4)*TEMP3)**2 320 CONTINUE ERK = ERK + (TEMP4*TEMP3)**2 330 CONTINUE IF (KM2 .EQ. 0) GO TO 340 C ......EXIT IF (KM2 .LT. 0) GO TO 350 ERKM2 = ABSH*SIG(KM1)*GSTR(KM2) 1 *DSQRT(ERKM2) 340 CONTINUE ERKM1 = ABSH*SIG(K)*GSTR(KM1)*DSQRT(ERKM1) 350 CONTINUE TEMP5 = ABSH*DSQRT(ERK) ERR = TEMP5*(G(K) - G(KP1)) ERK = TEMP5*SIG(KP1)*GSTR(K) KNEW = K C C TEST IF ORDER SHOULD BE LOWERED C IF (KM2 .EQ. 0) GO TO 370 C BEGIN BLOCK PERMITTING ...EXITS TO 360 C ...EXIT IF (KM2 .LT. 0) GO TO 360 IF (DMAX1(ERKM1,ERKM2) .LE. ERK) 1 KNEW = KM1 360 CONTINUE GO TO 380 370 CONTINUE IF (ERKM1 .LE. 0.5D0*ERK) KNEW = KM1 380 CONTINUE C C TEST IF STEP SUCCESSFUL C IF (ERR .GT. EPS) GO TO 520 C *** END BLOCK 3 *** C C *** BEGIN BLOCK 4 *** C THE STEP IS SUCCESSFUL. CORRECT THE C PREDICTED SOLUTION, EVALUATE THE DERIVATIVES C USING THE CORRECTED SOLUTION AND UPDATE THE C DIFFERENCES. DETERMINE BEST ORDER AND STEP C SIZE FOR NEXT STEP. C *** KOLD = K HOLD = H C C CORRECT AND EVALUATE C TEMP1 = H*G(KP1) IF (NORND) GO TO 400 DO 390 L = 1, NEQN RHO = TEMP1*(YP(L) - PHI(L,1)) 1 - PHI(L,16) Y(L) = P(L) + RHO PHI(L,15) = (Y(L) - P(L)) - RHO 390 CONTINUE GO TO 420 400 CONTINUE DO 410 L = 1, NEQN Y(L) = P(L) + TEMP1*(YP(L) - PHI(L,1)) 410 CONTINUE 420 CONTINUE CALL DF(X,Y,YP,RPAR,IPAR) C C UPDATE DIFFERENCES FOR NEXT STEP C DO 430 L = 1, NEQN PHI(L,KP1) = YP(L) - PHI(L,1) PHI(L,KP2) = PHI(L,KP1) - PHI(L,KP2) 430 CONTINUE DO 450 I = 1, K DO 440 L = 1, NEQN PHI(L,I) = PHI(L,I) + PHI(L,KP1) 440 CONTINUE 450 CONTINUE C C ESTIMATE ERROR AT ORDER K+1 UNLESS: C IN FIRST PHASE WHEN ALWAYS RAISE ORDER, C ALREADY DECIDED TO LOWER ORDER, C STEP SIZE NOT CONSTANT SO ESTIMATE C UNRELIABLE C ERKP1 = 0.0D0 IF (KNEW .EQ. KM1 .OR. K .EQ. 12) 1 PHASE1 = .FALSE. IF (PHASE1) GO TO 510 IF (KNEW .NE. KM1) GO TO 460 C C LOWER ORDER C K = KM1 ERK = ERKM1 C ..................EXIT GO TO 600 460 CONTINUE C ...............EXIT IF (KP1 .GT. NS) GO TO 600 DO 470 L = 1, NEQN ERKP1 = ERKP1 + (PHI(L,KP2)/WT(L))**2 470 CONTINUE ERKP1 = ABSH*GSTR(KP1)*DSQRT(ERKP1) C C USING ESTIMATED ERROR AT ORDER K+1, C DETERMINE APPROPRIATE ORDER FOR NEXT STEP C IF (K .GT. 1) GO TO 480 C ..................EXIT IF (ERKP1 .GE. 0.5D0*ERK) GO TO 600 GO TO 500 480 CONTINUE IF (ERKM1 .GT. DMIN1(ERK,ERKP1)) 1 GO TO 490 C C LOWER ORDER C K = KM1 ERK = ERKM1 C .....................EXIT GO TO 600 490 CONTINUE C ..................EXIT IF (ERKP1 .GE. ERK .OR. K .EQ. 12) 1 GO TO 600 500 CONTINUE 510 CONTINUE C C HERE ERKP1 .LT. ERK .LT. DMAX1(ERKM1,ERKM2) C ELSE ORDER WOULD HAVE BEEN LOWERED IN BLOCK C 2. THUS ORDER IS TO BE RAISED C C RAISE ORDER C K = KP1 ERK = ERKP1 C ............EXIT GO TO 600 520 CONTINUE C *** END BLOCK 2 *** C C *** BEGIN BLOCK 3 *** C THE STEP IS UNSUCCESSFUL. RESTORE X, PHI(*,*), C PSI(*) . IF THIRD CONSECUTIVE FAILURE, SET C ORDER TO ONE. IF STEP FAILS MORE THAN THREE C TIMES, CONSIDER AN OPTIMAL STEP SIZE. DOUBLE C ERROR TOLERANCE AND RETURN IF ESTIMATED STEP C SIZE IS TOO SMALL FOR MACHINE PRECISION. C *** C C RESTORE X, PHI(*,*) AND PSI(*) C PHASE1 = .FALSE. X = XOLD DO 540 I = 1, K TEMP1 = 1.0D0/BETA(I) IP1 = I + 1 DO 530 L = 1, NEQN PHI(L,I) = TEMP1*(PHI(L,I) - PHI(L,IP1)) 530 CONTINUE 540 CONTINUE IF (K .LT. 2) GO TO 560 DO 550 I = 2, K PSI(I-1) = PSI(I) - H 550 CONTINUE 560 CONTINUE C C ON THIRD FAILURE, SET ORDER TO ONE. THEREAFTER, C USE OPTIMAL STEP SIZE C IFAIL = IFAIL + 1 TEMP2 = 0.5D0 IF (IFAIL .NE. 3) GO TO 570 KNEW = 1 C ......EXIT GO TO 580 570 CONTINUE C ...EXIT IF (IFAIL .LT. 3) GO TO 580 IF (P5EPS .LT. 0.25D0*ERK) 1 TEMP2 = DSQRT(P5EPS/ERK) KNEW = 1 580 CONTINUE H = TEMP2*H K = KNEW IF (DABS(H) .GE. FOURU*DABS(X)) GO TO 590 CRASH = .TRUE. H = DSIGN(FOURU*DABS(X),H) EPS = EPS + EPS C ............EXIT GO TO 630 590 CONTINUE GO TO 70 600 CONTINUE C C WITH NEW ORDER DETERMINE APPROPRIATE STEP SIZE FOR NEXT C STEP C HNEW = H + H IF (PHASE1) GO TO 620 IF (P5EPS .GE. ERK*TWO(K+1)) GO TO 620 HNEW = H IF (P5EPS .GE. ERK) GO TO 610 TEMP2 = K + 1 R = (P5EPS/ERK)**(1.0D0/TEMP2) HNEW = ABSH*DMAX1(0.5D0,DMIN1(0.9D0,R)) HNEW = DSIGN(DMAX1(HNEW,FOURU*DABS(X)),H) 610 CONTINUE 620 CONTINUE H = HNEW 630 CONTINUE 640 CONTINUE 650 CONTINUE RETURN C *** END BLOCK 4 *** END SUBROUTINE DINTRP(X,Y,XOUT,YOUT,YPOUT,NEQN,KOLD,PHI,PSI) C***BEGIN PROLOGUE DINTRP C***DATE WRITTEN 820301 (YYMMDD) C***REVISION DATE 830701 (YYMMDD) C***CATEGORY NO. Z C***KEYWORDS ADAMS METHOD,DEPAC,INITIAL VALUE PROBLEMS,ODE, C ORDINARY DIFFERENTIAL EQUATIONS,PREDICTOR-CORRECTOR C***AUTHOR SHAMPINE L.F.(SNLA) C GORDON M.K. C***PURPOSE Approximates the solution at XOUT by evaluating the C polynomial computed in DSTEP2 at XOUT. Must be used C in conjunction with DSTEP2. C***DESCRIPTION C C Written by L. F. Shampine and M. K. Gordon C C Abstract C C C The methods in subroutine DSTEP2 approximate the solution near X C by a polynomial. Subroutine DINTRP approximates the solution at C XOUT by evaluating the polynomial there. Information defining this C polynomial is passed from DSTEP2 so DINTRP cannot be used alone. C C This code is completely explained and documented in the text, C COMPUTER SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS, THE INITIAL C VALUE PROBLEM by L. F. Shampine and M. K. Gordon. C Further details on use of this code are available in *SOLVING C ORDINARY DIFFERENTIAL EQUATIONS WITH ODE, STEP, AND DINTRP*, C by L. F. Shampine and M. K. Gordon, SLA-73-1060. C C Input to DINTRP -- C C The user provides storage in the calling program for the arrays in C the call list C DIMENSION Y(NEQN),YOUT(NEQN),YPOUT(NEQN),PHI(NEQN,16),PSI(12) C and defines C XOUT -- Point at which solution is desired. C The remaining parameters are defined in DSTEP2 and passed to C DINTRP from that subroutine C C Output from DINTRP -- C C YOUT(*) -- Solution at XOUT C YPOUT(*) -- Derivative of solution at XOUT C The remaining parameters are returned unaltered from their input C values. Integration with DSTEP2 may be continued. C***REFERENCES SHAMPINE, L. F., GORDON, M. K., *SOLVING ORDINARY C DIFFERENTIAL EQUATIONS WITH ODE, STEP, AND DINTRP*, C SLA-73-1060, SANDIA LABORATORIES, 1973. C***ROUTINES CALLED (NONE) C***END PROLOGUE DINTRP C INTEGER I, J, JM1, KI, KIP1, KOLD, L, LIMIT1, NEQN DOUBLE PRECISION ETA, G, GAMMA, HI, PHI, PSI, PSIJM1, RHO, TEMP1, 1 TEMP2, TEMP3, TERM, W, X, XOUT, Y, YOUT, YPOUT DIMENSION Y(NEQN),YOUT(NEQN),YPOUT(NEQN),PHI(NEQN,16),PSI(12) DIMENSION G(13),W(13),RHO(13) DATA G(1) /1.0D0/, RHO(1) /1.0D0/ C C***FIRST EXECUTABLE STATEMENT DINTRP HI = XOUT - X KI = KOLD + 1 KIP1 = KI + 1 C C INITIALIZE W(*) FOR COMPUTING G(*) C DO 10 I = 1, KI TEMP1 = I W(I) = 1.0D0/TEMP1 10 CONTINUE TERM = 0.0D0 C C COMPUTE G(*) C DO 30 J = 2, KI JM1 = J - 1 PSIJM1 = PSI(JM1) GAMMA = (HI + TERM)/PSIJM1 ETA = HI/PSIJM1 LIMIT1 = KIP1 - J DO 20 I = 1, LIMIT1 W(I) = GAMMA*W(I) - ETA*W(I+1) 20 CONTINUE G(J) = W(1) RHO(J) = GAMMA*RHO(JM1) TERM = PSIJM1 30 CONTINUE C C INTERPOLATE C DO 40 L = 1, NEQN YPOUT(L) = 0.0D0 YOUT(L) = 0.0D0 40 CONTINUE DO 60 J = 1, KI I = KIP1 - J TEMP2 = G(I) TEMP3 = RHO(I) DO 50 L = 1, NEQN YOUT(L) = YOUT(L) + TEMP2*PHI(L,I) YPOUT(L) = YPOUT(L) + TEMP3*PHI(L,I) 50 CONTINUE 60 CONTINUE DO 70 L = 1, NEQN YOUT(L) = Y(L) + HI*YOUT(L) 70 CONTINUE RETURN END SUBROUTINE DHSTRT(DF,NEQ,A,B,Y,YPRIME,ETOL,MORDER,SMALL,BIG,SPY, 1 PV,YP,SF,RPAR,IPAR,H) C***BEGIN PROLOGUE DHSTRT C***REFER TO DDEABM,DDEBDF,DDERKF C C DHSTRT computes a starting step size to be used in solving initial C value problems in ordinary differential equations. C C ********************************************************************** C ABSTRACT C C Subroutine DHSTRT computes a starting step size to be used by an C initial value method in solving ordinary differential equations. C It is based on an estimate of the local Lipschitz constant for the C differential equation (lower bound on a norm of the Jacobian) , C a bound on the differential equation (first derivative) , and C a bound on the partial derivative of the equation with respect to C the independent variable. C (all approximated near the initial point A) C C Subroutine DHSTRT uses a function subprogram DVNORM for computing C a vector norm. The maximum norm is presently utilized though it C can easily be replaced by any other vector norm. It is presumed C that any replacement norm routine would be carefully coded to C prevent unnecessary underflows or overflows from occurring, and C also, would not alter the vector or number of components. C C ********************************************************************** C On input you must provide the following C C DF -- This is a subroutine of the form C DF(X,U,UPRIME,RPAR,IPAR) C which defines the system of first order differential C equations to be solved. For the given values of X and the C vector U(*)=(U(1),U(2),...,U(NEQ)) , the subroutine must C evaluate the NEQ components of the system of differential C equations DU/DX=DF(X,U) and store the derivatives in the C array UPRIME(*), that is, UPRIME(I) = * DU(I)/DX * for C equations I=1,...,NEQ. C C Subroutine DF must not alter X or U(*). You must declare C the name DF in an external statement in your program that C calls DHSTRT. You must dimension U and UPRIME in DF. C C RPAR and IPAR are DOUBLE PRECISION and INTEGER parameter C arrays which you can use for communication between your C program and subroutine DF. They are not used or altered by C DHSTRT. If you do not need RPAR or IPAR, ignore these C parameters by treating them as dummy arguments. If you do C choose to use them, dimension them in your program and in C DF as arrays of appropriate length. C C NEQ -- This is the number of (first order) differential equations C to be integrated. C C A -- This is the initial point of integration. C C B -- This is a value of the independent variable used to define C the direction of integration. A reasonable choice is to C set B to the first point at which a solution is desired. C You can also use B, if necessary, to restrict the length C of the first integration step because the algorithm will C not compute a starting step length which is bigger than C DABS(B-A), unless B has been chosen too close to A. C (it is presumed that DHSTRT has been called with B C different from A on the machine being used. Also see the C discussion about the parameter SMALL.) C C Y(*) -- This is the vector of initial values of the NEQ solution C components at the initial point A. C C YPRIME(*) -- This is the vector of derivatives of the NEQ C solution components at the initial point A. C (defined by the differential equations in subroutine DF) C C ETOL -- This is the vector of error tolerances corresponding to C the NEQ solution components. It is assumed that all C elements are positive. Following the first integration C step, the tolerances are expected to be used by the C integrator in an error test which roughly requires that C DABS(LOCAL ERROR) .LE. ETOL C for each vector component. C C MORDER -- This is the order of the formula which will be used by C the initial value method for taking the first integration C step. C C SMALL -- This is a small positive machine dependent constant C which is used for protecting against computations with C numbers which are too small relative to the precision of C floating point arithmetic. SMALL should be set to C (approximately) the smallest positive DOUBLE PRECISION C number such that (1.+SMALL) .GT. 1. on the machine being C used. The quantity SMALL**(3/8) is used in computing C increments of variables for approximating derivatives by C differences. Also the algorithm will not compute a C starting step length which is smaller than C 100*SMALL*DABS(A). C C BIG -- This is a large positive machine dependent constant which C is used for preventing machine overflows. A reasonable C choice is to set big to (approximately) the square root of C the largest DOUBLE PRECISION number which can be held in C the machine. C C SPY(*),PV(*),YP(*),SF(*) -- These are DOUBLE PRECISION work C arrays of length NEQ which provide the routine with needed C storage space. C C RPAR,IPAR -- These are parameter arrays, of DOUBLE PRECISION and C INTEGER type, respectively, which can be used for C communication between your program and the DF subroutine. C They are not used or altered by DHSTRT. C C ********************************************************************** C On Output (after the return from DHSTRT), C C H -- is an appropriate starting step size to be attempted by the C differential equation method. C C All parameters in the call list remain unchanged except for C the working arrays SPY(*),PV(*),YP(*), and SF(*). C C ********************************************************************** C***ROUTINES CALLED DVNORM C***END PROLOGUE DHSTRT C INTEGER IPAR, J, K, LK, MIN0, MORDER, NEQ DOUBLE PRECISION A, ABSDX, B, BIG, DA, DABS, DBLE, DELF, DELY, 1 DFDUB, DFDXB, DLOG10, DMAX1, DMIN1, DSIGN, DSQRT, DVNORM, 2 DX, DY, ETOL, FBND, H, PV, RELPER, RPAR, SF, SMALL, SPY, 3 SRYDPB, TOLEXP, TOLMIN, TOLP, TOLSUM, Y, YDPB, YP, YPRIME DIMENSION Y(NEQ),YPRIME(NEQ),ETOL(NEQ),SPY(NEQ),PV(NEQ),YP(NEQ), 1 SF(NEQ),RPAR(1),IPAR(1) EXTERNAL DF C C .................................................................. C C BEGIN BLOCK PERMITTING ...EXITS TO 160 C***FIRST EXECUTABLE STATEMENT DHSTRT DX = B - A ABSDX = DABS(DX) RELPER = SMALL**0.375D0 C C ............................................................... C C COMPUTE AN APPROXIMATE BOUND (DFDXB) ON THE PARTIAL C DERIVATIVE OF THE EQUATION WITH RESPECT TO THE C INDEPENDENT VARIABLE. PROTECT AGAINST AN OVERFLOW. C ALSO COMPUTE A BOUND (FBND) ON THE FIRST DERIVATIVE C LOCALLY. C DA = DSIGN(DMAX1(DMIN1(RELPER*DABS(A),ABSDX), 1 100.0D0*SMALL*DABS(A)),DX) IF (DA .EQ. 0.0D0) DA = RELPER*DX CALL DF(A+DA,Y,SF,RPAR,IPAR) DO 10 J = 1, NEQ YP(J) = SF(J) - YPRIME(J) 10 CONTINUE DELF = DVNORM(YP,NEQ) DFDXB = BIG IF (DELF .LT. BIG*DABS(DA)) DFDXB = DELF/DABS(DA) FBND = DVNORM(SF,NEQ) C C ............................................................... C C COMPUTE AN ESTIMATE (DFDUB) OF THE LOCAL LIPSCHITZ C CONSTANT FOR THE SYSTEM OF DIFFERENTIAL EQUATIONS. THIS C ALSO REPRESENTS AN ESTIMATE OF THE NORM OF THE JACOBIAN C LOCALLY. THREE ITERATIONS (TWO WHEN NEQ=1) ARE USED TO C ESTIMATE THE LIPSCHITZ CONSTANT BY NUMERICAL DIFFERENCES. C THE FIRST PERTURBATION VECTOR IS BASED ON THE INITIAL C DERIVATIVES AND DIRECTION OF INTEGRATION. THE SECOND C PERTURBATION VECTOR IS FORMED USING ANOTHER EVALUATION OF C THE DIFFERENTIAL EQUATION. THE THIRD PERTURBATION VECTOR C IS FORMED USING PERTURBATIONS BASED ONLY ON THE INITIAL C VALUES. COMPONENTS THAT ARE ZERO ARE ALWAYS CHANGED TO C NON-ZERO VALUES (EXCEPT ON THE FIRST ITERATION). WHEN C INFORMATION IS AVAILABLE, CARE IS TAKEN TO ENSURE THAT C COMPONENTS OF THE PERTURBATION VECTOR HAVE SIGNS WHICH ARE C CONSISTENT WITH THE SLOPES OF LOCAL SOLUTION CURVES. C ALSO CHOOSE THE LARGEST BOUND (FBND) FOR THE FIRST C DERIVATIVE. C C PERTURBATION VECTOR SIZE IS HELD C CONSTANT FOR ALL ITERATIONS. COMPUTE C THIS CHANGE FROM THE C SIZE OF THE VECTOR OF INITIAL C VALUES. DELY = RELPER*DVNORM(Y,NEQ) IF (DELY .EQ. 0.0D0) DELY = RELPER DELY = DSIGN(DELY,DX) DELF = DVNORM(YPRIME,NEQ) FBND = DMAX1(FBND,DELF) IF (DELF .EQ. 0.0D0) GO TO 30 C USE INITIAL DERIVATIVES FOR FIRST PERTURBATION DO 20 J = 1, NEQ SPY(J) = YPRIME(J) YP(J) = YPRIME(J) 20 CONTINUE GO TO 50 30 CONTINUE C CANNOT HAVE A NULL PERTURBATION VECTOR DO 40 J = 1, NEQ SPY(J) = 0.0D0 YP(J) = 1.0D0 40 CONTINUE DELF = DVNORM(YP,NEQ) 50 CONTINUE C DFDUB = 0.0D0 LK = MIN0(NEQ+1,3) DO 140 K = 1, LK C DEFINE PERTURBED VECTOR OF INITIAL VALUES DO 60 J = 1, NEQ PV(J) = Y(J) + DELY*(YP(J)/DELF) 60 CONTINUE IF (K .EQ. 2) GO TO 80 C EVALUATE DERIVATIVES ASSOCIATED WITH PERTURBED C VECTOR AND COMPUTE CORRESPONDING DIFFERENCES CALL DF(A,PV,YP,RPAR,IPAR) DO 70 J = 1, NEQ PV(J) = YP(J) - YPRIME(J) 70 CONTINUE GO TO 100 80 CONTINUE C USE A SHIFTED VALUE OF THE INDEPENDENT VARIABLE C IN COMPUTING ONE ESTIMATE CALL DF(A+DA,PV,YP,RPAR,IPAR) DO 90 J = 1, NEQ PV(J) = YP(J) - SF(J) 90 CONTINUE 100 CONTINUE C CHOOSE LARGEST BOUNDS ON THE FIRST DERIVATIVE C AND A LOCAL LIPSCHITZ CONSTANT FBND = DMAX1(FBND,DVNORM(YP,NEQ)) DELF = DVNORM(PV,NEQ) C ...EXIT IF (DELF .GE. BIG*DABS(DELY)) GO TO 150 DFDUB = DMAX1(DFDUB,DELF/DABS(DELY)) C ......EXIT IF (K .EQ. LK) GO TO 160 C CHOOSE NEXT PERTURBATION VECTOR IF (DELF .EQ. 0.0D0) DELF = 1.0D0 DO 130 J = 1, NEQ IF (K .EQ. 2) GO TO 110 DY = DABS(PV(J)) IF (DY .EQ. 0.0D0) DY = DELF GO TO 120 110 CONTINUE DY = Y(J) IF (DY .EQ. 0.0D0) DY = DELY/RELPER 120 CONTINUE IF (SPY(J) .EQ. 0.0D0) SPY(J) = YP(J) IF (SPY(J) .NE. 0.0D0) DY = DSIGN(DY,SPY(J)) YP(J) = DY 130 CONTINUE DELF = DVNORM(YP,NEQ) 140 CONTINUE 150 CONTINUE C C PROTECT AGAINST AN OVERFLOW DFDUB = BIG 160 CONTINUE C C .................................................................. C C COMPUTE A BOUND (YDPB) ON THE NORM OF THE SECOND DERIVATIVE C YDPB = DFDXB + DFDUB*FBND C C .................................................................. C C DEFINE THE TOLERANCE PARAMETER UPON WHICH THE STARTING STEP C SIZE IS TO BE BASED. A VALUE IN THE MIDDLE OF THE ERROR C TOLERANCE RANGE IS SELECTED. C TOLMIN = BIG TOLSUM = 0.0D0 DO 170 K = 1, NEQ TOLEXP = DLOG10(ETOL(K)) TOLMIN = DMIN1(TOLMIN,TOLEXP) TOLSUM = TOLSUM + TOLEXP 170 CONTINUE TOLP = 10.0D0 1 **(0.5D0*(TOLSUM/DBLE(FLOAT(NEQ)) + TOLMIN) 2 /DBLE(FLOAT(MORDER+1))) C C .................................................................. C C COMPUTE A STARTING STEP SIZE BASED ON THE ABOVE FIRST AND C SECOND DERIVATIVE INFORMATION C C RESTRICT THE STEP LENGTH TO BE NOT BIGGER C THAN DABS(B-A). (UNLESS B IS TOO CLOSE C TO A) H = ABSDX C IF (YDPB .NE. 0.0D0 .OR. FBND .NE. 0.0D0) GO TO 180 C C BOTH FIRST DERIVATIVE TERM (FBND) AND SECOND C DERIVATIVE TERM (YDPB) ARE ZERO IF (TOLP .LT. 1.0D0) H = ABSDX*TOLP GO TO 200 180 CONTINUE C IF (YDPB .NE. 0.0D0) GO TO 190 C C ONLY SECOND DERIVATIVE TERM (YDPB) IS ZERO IF (TOLP .LT. FBND*ABSDX) H = TOLP/FBND GO TO 200 190 CONTINUE C C SECOND DERIVATIVE TERM (YDPB) IS NON-ZERO SRYDPB = DSQRT(0.5D0*YDPB) IF (TOLP .LT. SRYDPB*ABSDX) H = TOLP/SRYDPB 200 CONTINUE C C FURTHER RESTRICT THE STEP LENGTH TO BE NOT C BIGGER THAN 1/DFDUB IF (H*DFDUB .GT. 1.0D0) H = 1.0D0/DFDUB C C FINALLY, RESTRICT THE STEP LENGTH TO BE NOT C SMALLER THAN 100*SMALL*DABS(A). HOWEVER, IF C A=0. AND THE COMPUTED H UNDERFLOWED TO ZERO, C THE ALGORITHM RETURNS SMALL*DABS(B) FOR THE C STEP LENGTH. H = DMAX1(H,100.0D0*SMALL*DABS(A)) IF (H .EQ. 0.0D0) H = SMALL*DABS(B) C C NOW SET DIRECTION OF INTEGRATION H = DSIGN(H,DX) C RETURN END DOUBLE PRECISION FUNCTION DVNORM(V,NCOMP) C***BEGIN PROLOGUE DVNORM C***REFER TO DDEABM,DDEBDF,DDERKF C C Compute the maximum norm of the vector V(*) of length NCOMP and C return the result as DVNORM C***ROUTINES CALLED (NONE) C***END PROLOGUE DVNORM C INTEGER K, NCOMP DOUBLE PRECISION DABS, DMAX1, V DIMENSION V(NCOMP) C***FIRST EXECUTABLE STATEMENT DVNORM DVNORM = 0.0D0 DO 10 K = 1, NCOMP DVNORM = DMAX1(DVNORM,DABS(V(K))) 10 CONTINUE RETURN END DOUBLE PRECISION FUNCTION D1MACH(I) C***BEGIN PROLOGUE D1MACH C***DATE WRITTEN 750101 (YYMMDD) C***REVISION DATE 831014 (YYMMDD) C***CATEGORY NO. R1 C***KEYWORDS MACHINE CONSTANTS C***AUTHOR FOX, P. A., (BELL LABS) C HALL, A. D., (BELL LABS) C SCHRYER, N. L., (BELL LABS) C***PURPOSE Returns double precision machine dependent constants C***DESCRIPTION C C D1MACH can be used to obtain machine-dependent parameters C for the local machine environment. It is a function C subprogram with one (input) argument, and can be called C as follows, for example C C D = D1MACH(I) C C where I=1,...,5. The (output) value of D above is C determined by the (input) value of I. The results for C various values of I are discussed below. C C Double-precision machine constants C D1MACH( 1) = B**(EMIN-1), the smallest positive magnitude. C D1MACH( 2) = B**EMAX*(1 - B**(-T)), the largest magnitude. C D1MACH( 3) = B**(-T), the smallest relative spacing. C D1MACH( 4) = B**(1-T), the largest relative spacing. C D1MACH( 5) = LOG10(B) C***REFERENCES FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A C PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188. C***END PROLOGUE D1MACH C INTEGER SMALL(4) INTEGER LARGE(4) INTEGER RIGHT(4) INTEGER DIVER(4) INTEGER LOG10(4) C DOUBLE PRECISION DMACH(5) C EQUIVALENCE (DMACH(1),SMALL(1)) EQUIVALENCE (DMACH(2),LARGE(1)) EQUIVALENCE (DMACH(3),RIGHT(1)) EQUIVALENCE (DMACH(4),DIVER(1)) EQUIVALENCE (DMACH(5),LOG10(1)) C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C DATA SMALL(1) / ZC00800000 / C DATA SMALL(2) / Z000000000 / C C DATA LARGE(1) / ZDFFFFFFFF / C DATA LARGE(2) / ZFFFFFFFFF / C C DATA RIGHT(1) / ZCC5800000 / C DATA RIGHT(2) / Z000000000 / C C DATA DIVER(1) / ZCC6800000 / C DATA DIVER(2) / Z000000000 / C C DATA LOG10(1) / ZD00E730E7 / C DATA LOG10(2) / ZC77800DC0 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. C C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O0000000000000000 / C C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O0007777777777777 / C C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 / C C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. C C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O7770000000000000 / C C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O7777777777777777 / C C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES. C C DATA SMALL(1) / 00564000000000000000B / C DATA SMALL(2) / 00000000000000000000B / C C DATA LARGE(1) / 37757777777777777777B / C DATA LARGE(2) / 37157777777777777777B / C C DATA RIGHT(1) / 15624000000000000000B / C DATA RIGHT(2) / 00000000000000000000B / C C DATA DIVER(1) / 15634000000000000000B / C DATA DIVER(2) / 00000000000000000000B / C C DATA LOG10(1) / 17164642023241175717B / C DATA LOG10(2) / 16367571421742254654B / C C MACHINE CONSTANTS FOR THE CRAY 1 C C DATA SMALL(1) / 201354000000000000000B / C DATA SMALL(2) / 000000000000000000000B / C C DATA LARGE(1) / 577767777777777777777B / C DATA LARGE(2) / 000007777777777777774B / C C DATA RIGHT(1) / 376434000000000000000B / C DATA RIGHT(2) / 000000000000000000000B / C C DATA DIVER(1) / 376444000000000000000B / C DATA DIVER(2) / 000000000000000000000B / C C DATA LOG10(1) / 377774642023241175717B / C DATA LOG10(2) / 000007571421742254654B / C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD - C STATIC DMACH(5) C C DATA SMALL/20K,3*0/,LARGE/77777K,3*177777K/ C DATA RIGHT/31420K,3*0/,DIVER/32020K,3*0/ C DATA LOG10/40423K,42023K,50237K,74776K/ C C MACHINE CONSTANTS FOR THE HARRIS 220 C C DATA SMALL(1),SMALL(2) / '20000000, '00000201 / C DATA LARGE(1),LARGE(2) / '37777777, '37777577 / C DATA RIGHT(1),RIGHT(2) / '20000000, '00000333 / C DATA DIVER(1),DIVER(2) / '20000000, '00000334 / C DATA LOG10(1),LOG10(2) / '23210115, '10237777 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES. C C DATA SMALL(1),SMALL(2) / O402400000000, O000000000000 / C DATA LARGE(1),LARGE(2) / O376777777777, O777777777777 / C DATA RIGHT(1),RIGHT(2) / O604400000000, O000000000000 / C DATA DIVER(1),DIVER(2) / O606400000000, O000000000000 / C DATA LOG10(1),LOG10(2) / O776464202324, O117571775714 / C C MACHINE CONSTANTS FOR THE HP 2100 C THREE WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA SMALL(1), SMALL(2), SMALL(3) / 40000B, 0, 1 / C DATA LARGE(1), LARGE(2), LARGE(3) / 77777B, 177777B, 177776B / C DATA RIGHT(1), RIGHT(2), RIGHT(3) / 40000B, 0, 265B / C DATA DIVER(1), DIVER(2), DIVER(3) / 40000B, 0, 276B / C DATA LOG10(1), LOG10(2), LOG10(3) / 46420B, 46502B, 77777B / C C C MACHINE CONSTANTS FOR THE HP 2100 C FOUR WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA SMALL(1), SMALL(2) / 40000B, 0 / C DATA SMALL(3), SMALL(4) / 0, 1 / C DATA LARGE(1), LARGE(2) / 77777B, 177777B / C DATA LARGE(3), LARGE(4) / 177777B, 177776B / C DATA RIGHT(1), RIGHT(2) / 40000B, 0 / C DATA RIGHT(3), RIGHT(4) / 0, 225B / C DATA DIVER(1), DIVER(2) / 40000B, 0 / C DATA DIVER(3), DIVER(4) / 0, 227B / C DATA LOG10(1), LOG10(2) / 46420B, 46502B / C DATA LOG10(3), LOG10(4) / 76747B, 176377B / C C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND C THE PERKIN ELMER (INTERDATA) 7/32. C C DATA SMALL(1),SMALL(2) / Z00100000, Z00000000 / C DATA LARGE(1),LARGE(2) / Z7FFFFFFF, ZFFFFFFFF / C DATA RIGHT(1),RIGHT(2) / Z33100000, Z00000000 / C DATA DIVER(1),DIVER(2) / Z34100000, Z00000000 / C DATA LOG10(1),LOG10(2) / Z41134413, Z509F79FF / C C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR). C C DATA SMALL(1),SMALL(2) / "033400000000, "000000000000 / C DATA LARGE(1),LARGE(2) / "377777777777, "344777777777 / C DATA RIGHT(1),RIGHT(2) / "113400000000, "000000000000 / C DATA DIVER(1),DIVER(2) / "114400000000, "000000000000 / C DATA LOG10(1),LOG10(2) / "177464202324, "144117571776 / C C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR). C C DATA SMALL(1),SMALL(2) / "000400000000, "000000000000 / C DATA LARGE(1),LARGE(2) / "377777777777, "377777777777 / C DATA RIGHT(1),RIGHT(2) / "103400000000, "000000000000 / C DATA DIVER(1),DIVER(2) / "104400000000, "000000000000 / C DATA LOG10(1),LOG10(2) / "177464202324, "476747767461 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1),SMALL(2) / 8388608, 0 / C DATA LARGE(1),LARGE(2) / 2147483647, -1 / C DATA RIGHT(1),RIGHT(2) / 612368384, 0 / C DATA DIVER(1),DIVER(2) / 620756992, 0 / C DATA LOG10(1),LOG10(2) / 1067065498, -2063872008 / C C DATA SMALL(1),SMALL(2) / O00040000000, O00000000000 / C DATA LARGE(1),LARGE(2) / O17777777777, O37777777777 / C DATA RIGHT(1),RIGHT(2) / O04440000000, O00000000000 / C DATA DIVER(1),DIVER(2) / O04500000000, O00000000000 / C DATA LOG10(1),LOG10(2) / O07746420232, O20476747770 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1),SMALL(2) / 128, 0 / C DATA SMALL(3),SMALL(4) / 0, 0 / C C DATA LARGE(1),LARGE(2) / 32767, -1 / C DATA LARGE(3),LARGE(4) / -1, -1 / C C DATA RIGHT(1),RIGHT(2) / 9344, 0 / C DATA RIGHT(3),RIGHT(4) / 0, 0 / C C DATA DIVER(1),DIVER(2) / 9472, 0 / C DATA DIVER(3),DIVER(4) / 0, 0 / C C DATA LOG10(1),LOG10(2) / 16282, 8346 / C DATA LOG10(3),LOG10(4) / -31493, -12296 / C C DATA SMALL(1),SMALL(2) / O000200, O000000 / C DATA SMALL(3),SMALL(4) / O000000, O000000 / C C DATA LARGE(1),LARGE(2) / O077777, O177777 / C DATA LARGE(3),LARGE(4) / O177777, O177777 / C C DATA RIGHT(1),RIGHT(2) / O022200, O000000 / C DATA RIGHT(3),RIGHT(4) / O000000, O000000 / C C DATA DIVER(1),DIVER(2) / O022400, O000000 / C DATA DIVER(3),DIVER(4) / O000000, O000000 / C C DATA LOG10(1),LOG10(2) / O037632, O020232 / C DATA LOG10(3),LOG10(4) / O102373, O147770 / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILER C C DATA SMALL(1),SMALL(2) / O000040000000, O000000000000 / C DATA LARGE(1),LARGE(2) / O377777777777, O777777777777 / C DATA RIGHT(1),RIGHT(2) / O170540000000, O000000000000 / C DATA DIVER(1),DIVER(2) / O170640000000, O000000000000 / C DATA LOG10(1),LOG10(2) / O177746420232, O411757177572 / C C C MACHINE CONSTANTS FOR VAX 11/780 C (EXPRESSED IN INTEGER AND HEXADECIMAL) C ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS*** C *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS*** C C DATA SMALL(1), SMALL(2) / 128, 0 / C DATA LARGE(1), LARGE(2) / -32769, -1 / C DATA RIGHT(1), RIGHT(2) / 9344, 0 / C DATA DIVER(1), DIVER(2) / 9472, 0 / C DATA LOG10(1), LOG10(2) / 546979738, -805796613 / C C DATA SMALL(1), SMALL(2) / Z00000080, Z00000000 / C DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / C DATA RIGHT(1), RIGHT(2) / Z00002480, Z00000000 / C DATA DIVER(1), DIVER(2) / Z00002500, Z00000000 / C DATA LOG10(1), LOG10(2) / Z209A3F9A, ZCFF884FB / C C MACHINE CONSTANTS FOR VAX 11/780 (G-FLOATING) C (EXPRESSED IN INTEGER AND HEXADECIMAL) C ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS*** C *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS*** C C DATA SMALL(1), SMALL(2) / 16, 0 / C DATA LARGE(1), LARGE(2) / -32769, -1 / C DATA RIGHT(1), RIGHT(2) / 15552, 0 / C DATA DIVER(1), DIVER(2) / 15568, 0 / C DATA LOG10(1), LOG10(2) / 1142112243, 2046775455 / C C DATA SMALL(1), SMALL(2) / Z00000010, Z00000000 / C DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / C DATA RIGHT(1), RIGHT(2) / Z00003CC0, Z00000000 / C DATA DIVER(1), DIVER(2) / Z00003CD0, Z00000000 / C DATA LOG10(1), LOG10(2) / Z44133FF3, Z79FF509F / C C MACHINE CONSTANTS FOR IBM/AT OR COMPATIBILES AND C Microsoft (R) FORTRAN Optimizing Compiler Version 4.00 C Copyright (c) Microsoft Corp 1987. All rights reserved. C (ADDED BY D.KEKEZ 88/10/04 YY/MM/DD) C DATA SMALL(1), SMALL(2) / 0 , 1048576/ DATA LARGE(1), LARGE(2) /-1 , 2147483647/ DATA RIGHT(1), RIGHT(2) / 0 , 1017118720/ DATA DIVER(1), DIVER(2) / 0 , 1018167296/ DATA LOG10(1), LOG10(2) / 1352628735, 1070810131/ C C C***FIRST EXECUTABLE STATEMENT D1MACH IF (I .LT. 1 .OR. I .GT. 5)THEN WRITE(6,1000)I 1000 FORMAT(' D1MACH - ARGUMENT I IS OUT OF BOUNDS, I=',I6) STOP ENDIF D1MACH = DMACH(I) RETURN C END SUBROUTINE XERRWV(MESSG,NMESSG,NERR,LEVEL,NI,I1,I2,NR,R1,R2) C***BEGIN PROLOGUE XERRWV C***DATE WRITTEN 800319 (YYMMDD) C***REVISION DATE 840405 (YYMMDD) C***CATEGORY NO. R3C C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Processes error message allowing 2 integer and two real C values to be included in the message. C***DESCRIPTION C C Abstract C XERRWV processes a diagnostic message, in a manner C determined by the value of LEVEL and the current value C of the library error control flag, KONTRL. C (See subroutine XSETF for details.) C In addition, up to two integer values and two real C values may be printed along with the message. C C Description of Parameters C --Input-- C MESSG - the Hollerith message to be processed. C NMESSG- the actual number of characters in MESSG. C NERR - the error number associated with this message. C NERR must not be zero. C LEVEL - error category. C =2 means this is an unconditionally fatal error. C =1 means this is a recoverable error. (I.e., it is C non-fatal if XSETF has been appropriately called.) C =0 means this is a warning message only. C =-1 means this is a warning message which is to be C printed at most once, regardless of how many C times this call is executed. C NI - number of integer values to be printed. (0 to 2) C I1 - first integer value. C I2 - second integer value. C NR - number of real values to be printed. (0 to 2) C R1 - first real value. C R2 - second real value. C C Examples C CALL XERROR(29HSMOOTH -- NUM (=I1) WAS ZERO.,29,1,2, C 1 1,NUM,0,0,0.,0.) C CALL XERRWV(54HQUADXY -- REQUESTED ERROR (R1) LESS THAN MINIMUM C 1 (R2).,54,77,1,0,0,0,2,ERRREQ,ERRMIN) C C Written by Ron Jones, with SLATEC Common Math Library Subcommittee C Latest revision --- 19 Mar 1980 C***REFERENCES JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR C HANDLING PACKAGE*, SAND78-1189, SANDIA LABORATORIES, C 1978. C***ROUTINES CALLED FDUMP,I1MACH,J4SAVE,XERABT,XERCTL,XERPRT,XERSAV, C XGETUA C***END PROLOGUE XERRWV C DIMENSION MESSG(NMESSG),LUN(5) C***FIRST EXECUTABLE STATEMENT XERRWV LKNTRL = J4SAVE(2,0,.FALSE.) MAXMES = J4SAVE(4,0,.FALSE.) C CHECK FOR VALID INPUT IF ((NMESSG.GT.0).AND.(NERR.NE.0).AND. 1 (LEVEL.GE.(-1)).AND.(LEVEL.LE.2)) GO TO 10 IF (LKNTRL.GT.0) CALL XERPRT(17HFATAL ERROR IN...,17) CALL XERPRT(23HXERROR -- INVALID INPUT,23) IF (LKNTRL.GT.0) CALL FDUMP IF (LKNTRL.GT.0) CALL XERPRT(29HJOB ABORT DUE TO FATAL ERROR., 1 29) IF (LKNTRL.GT.0) CALL XERSAV(1H ,0,0,0,KDUMMY) CALL XERABT(23HXERROR -- INVALID INPUT,23) RETURN 10 CONTINUE C RECORD MESSAGE JUNK = J4SAVE(1,NERR,.TRUE.) CALL XERSAV(MESSG,NMESSG,NERR,LEVEL,KOUNT) C LET USER OVERRIDE LFIRST = MESSG(1) LMESSG = NMESSG LERR = NERR LLEVEL = LEVEL CALL XERCTL(LFIRST,LMESSG,LERR,LLEVEL,LKNTRL) C RESET TO ORIGINAL VALUES LMESSG = NMESSG LERR = NERR LLEVEL = LEVEL LKNTRL = MAX0(-2,MIN0(2,LKNTRL)) MKNTRL = IABS(LKNTRL) C DECIDE WHETHER TO PRINT MESSAGE IF ((LLEVEL.LT.2).AND.(LKNTRL.EQ.0)) GO TO 100 IF (((LLEVEL.EQ.(-1)).AND.(KOUNT.GT.MIN0(1,MAXMES))) 1.OR.((LLEVEL.EQ.0) .AND.(KOUNT.GT.MAXMES)) 2.OR.((LLEVEL.EQ.1) .AND.(KOUNT.GT.MAXMES).AND.(MKNTRL.EQ.1)) 3.OR.((LLEVEL.EQ.2) .AND.(KOUNT.GT.MAX0(1,MAXMES)))) GO TO 100 IF (LKNTRL.LE.0) GO TO 20 CALL XERPRT(1H ,1) C INTRODUCTION IF (LLEVEL.EQ.(-1)) CALL XERPRT 1(57HWARNING MESSAGE...THIS MESSAGE WILL ONLY BE PRINTED ONCE.,57) IF (LLEVEL.EQ.0) CALL XERPRT(13HWARNING IN...,13) IF (LLEVEL.EQ.1) CALL XERPRT 1 (23HRECOVERABLE ERROR IN...,23) IF (LLEVEL.EQ.2) CALL XERPRT(17HFATAL ERROR IN...,17) 20 CONTINUE C MESSAGE CALL XERPRT(MESSG,LMESSG) CALL XGETUA(LUN,NUNIT) DO 50 KUNIT=1,NUNIT IUNIT = LUN(KUNIT) IF (IUNIT.EQ.0) IUNIT = I1MACH(4) IF (NI.GE.1) WRITE (IUNIT,22) I1 IF (NI.GE.2) WRITE (IUNIT,23) I2 IF (NR.GE.1) WRITE (IUNIT,24) R1 IF (NR.GE.2) WRITE (IUNIT,25) R2 22 FORMAT (11X, 'IN ABOVE MESSAGE, I1=',I10) 23 FORMAT (11X, 'IN ABOVE MESSAGE, I2=',I10) 24 FORMAT (11X, 'IN ABOVE MESSAGE, R1=',E20.10) 25 FORMAT (11X, 'IN ABOVE MESSAGE, R2=',E20.10) IF (LKNTRL.LE.0) GO TO 40 C ERROR NUMBER WRITE (IUNIT,30) LERR 30 FORMAT (' ERROR NUMBER =',I10) 40 CONTINUE 50 CONTINUE C TRACE-BACK CALL FDUMP 100 CONTINUE IFATAL = 0 IF ((LLEVEL.EQ.2).OR.((LLEVEL.EQ.1).AND.(MKNTRL.EQ.2))) 1IFATAL = 1 C QUIT HERE IF MESSAGE IS NOT FATAL IF (IFATAL.LE.0) RETURN IF ((LKNTRL.LE.0).OR.(KOUNT.GT.MAX0(1,MAXMES))) GO TO 120 C PRINT REASON FOR ABORT IF (LLEVEL.EQ.1) CALL XERPRT 1 (35HJOB ABORT DUE TO UNRECOVERED ERROR.,35) IF (LLEVEL.EQ.2) CALL XERPRT 1 (29HJOB ABORT DUE TO FATAL ERROR.,29) C PRINT ERROR SUMMARY CALL XERSAV(1H ,-1,0,0,KDUMMY) 120 CONTINUE C ABORT IF ((LLEVEL.EQ.2).AND.(KOUNT.GT.MAX0(1,MAXMES))) LMESSG = 0 CALL XERABT(MESSG,LMESSG) RETURN END SUBROUTINE FDUMP C***BEGIN PROLOGUE FDUMP C***DATE WRITTEN 790801 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. Z C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Symbolic dump (should be locally written) C***DESCRIPTION C C ***Note*** Machine Dependent Routine C FDUMP is intended to be replaced by a locally written C version which produces a symbolic dump. Failing this, C it should be replaced by a version which prints the C subprogram nesting list. Note that this dump must be C printed on each of up to five files, as indicated by the C XGETUA routine. See XSETUA and XGETUA for details. C C Written by Ron Jones, with SLATEC Common Math Library Subcommittee C Latest revision --- 23 May 1979 C***REFERENCES JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR C HANDLING PACKAGE*, SAND78-1189, SANDIA LABORATORIES, C 1978. C***ROUTINES CALLED (NONE) C***END PROLOGUE FDUMP C C***FIRST EXECUTABLE STATEMENT FDUMP RETURN END INTEGER FUNCTION I1MACH(I) C***BEGIN PROLOGUE I1MACH C***DATE WRITTEN 750101 (YYMMDD) C***REVISION DATE 840405 (YYMMDD) C***CATEGORY NO. R1 C***KEYWORDS MACHINE CONSTANTS C***AUTHOR FOX, P. A., (BELL LABS) C HALL, A. D., (BELL LABS) C SCHRYER, N. L., (BELL LABS) C***PURPOSE Returns integer machine dependent constants C***DESCRIPTION C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C These machine constant routines must be activated for C a particular environment. C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C I1MACH can be used to obtain machine-dependent parameters C for the local machine environment. It is a function C subroutine with one (input) argument, and can be called C as follows, for example C C K = I1MACH(I) C C where I=1,...,16. The (output) value of K above is C determined by the (input) value of I. The results for C various values of I are discussed below. C C I/O unit numbers. C I1MACH( 1) = the standard input unit. C I1MACH( 2) = the standard output unit. C I1MACH( 3) = the standard punch unit. C I1MACH( 4) = the standard error message unit. C C Words. C I1MACH( 5) = the number of bits per integer storage unit. C I1MACH( 6) = the number of characters per integer storage unit. C C Integers. C assume integers are represented in the S-digit, base-A form C C sign ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) ) C C where 0 .LE. X(I) .LT. A for I=0,...,S-1. C I1MACH( 7) = A, the base. C I1MACH( 8) = S, the number of base-A digits. C I1MACH( 9) = A**S - 1, the largest magnitude. C C Floating-Point Numbers. C Assume floating-point numbers are represented in the T-digit, C base-B form C sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) ) C C where 0 .LE. X(I) .LT. B for I=1,...,T, C 0 .LT. X(1), and EMIN .LE. E .LE. EMAX. C I1MACH(10) = B, the base. C C Single-Precision C I1MACH(11) = T, the number of base-B digits. C I1MACH(12) = EMIN, the smallest exponent E. C I1MACH(13) = EMAX, the largest exponent E. C C Double-Precision C I1MACH(14) = T, the number of base-B digits. C I1MACH(15) = EMIN, the smallest exponent E. C I1MACH(16) = EMAX, the largest exponent E. C C To alter this function for a particular environment, C the desired set of DATA statements should be activated by C removing the C from column 1. Also, the values of C I1MACH(1) - I1MACH(4) should be checked for consistency C with the local operating system. C***REFERENCES FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A C PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188. C***ROUTINES CALLED (NONE) C***END PROLOGUE I1MACH C INTEGER IMACH(16),OUTPUT EQUIVALENCE (IMACH(4),OUTPUT) C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C DATA IMACH( 1) / 7 / C DATA IMACH( 2) / 2 / C DATA IMACH( 3) / 2 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 33 / C DATA IMACH( 9) / Z1FFFFFFFF / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -256 / C DATA IMACH(13) / 255 / C DATA IMACH(14) / 60 / C DATA IMACH(15) / -256 / C DATA IMACH(16) / 255 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / 6 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 39 / C DATA IMACH( 9) / O0007777777777777 / C DATA IMACH(10) / 8 / C DATA IMACH(11) / 13 / C DATA IMACH(12) / -50 / C DATA IMACH(13) / 76 / C DATA IMACH(14) / 26 / C DATA IMACH(15) / -50 / C DATA IMACH(16) / 76 / C C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / 6 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 39 / C DATA IMACH( 9) / O0007777777777777 / C DATA IMACH(10) / 8 / C DATA IMACH(11) / 13 / C DATA IMACH(12) / -50 / C DATA IMACH(13) / 76 / C DATA IMACH(14) / 26 / C DATA IMACH(15) / -32754 / C DATA IMACH(16) / 32780 / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) /6LOUTPUT/ C DATA IMACH( 5) / 60 / C DATA IMACH( 6) / 10 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 48 / C DATA IMACH( 9) / 00007777777777777777B / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 47 / C DATA IMACH(12) / -929 / C DATA IMACH(13) / 1070 / C DATA IMACH(14) / 94 / C DATA IMACH(15) / -929 / C DATA IMACH(16) / 1069 / C C MACHINE CONSTANTS FOR THE CRAY 1 C C DATA IMACH( 1) / 100 / C DATA IMACH( 2) / 101 / C DATA IMACH( 3) / 102 / C DATA IMACH( 4) / 101 / C DATA IMACH( 5) / 64 / C DATA IMACH( 6) / 8 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 63 / C DATA IMACH( 9) / 777777777777777777777B / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 47 / C DATA IMACH(12) / -8189 / C DATA IMACH(13) / 8190 / C DATA IMACH(14) / 94 / C DATA IMACH(15) / -8099 / C DATA IMACH(16) / 8190 / C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C DATA IMACH( 1) / 11 / C DATA IMACH( 2) / 12 / C DATA IMACH( 3) / 8 / C DATA IMACH( 4) / 10 / C DATA IMACH( 5) / 16 / C DATA IMACH( 6) / 2 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 15 / C DATA IMACH( 9) /32767 / C DATA IMACH(10) / 16 / C DATA IMACH(11) / 6 / C DATA IMACH(12) / -64 / C DATA IMACH(13) / 63 / C DATA IMACH(14) / 14 / C DATA IMACH(15) / -64 / C DATA IMACH(16) / 63 / C C MACHINE CONSTANTS FOR THE HARRIS 220 C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 0 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / 3 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 23 / C DATA IMACH( 9) / 8388607 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 23 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 38 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 43 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 6 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / O377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 63 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE HP 2100 C 3 WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA IMACH(1) / 5/ C DATA IMACH(2) / 6 / C DATA IMACH(3) / 4 / C DATA IMACH(4) / 1 / C DATA IMACH(5) / 16 / C DATA IMACH(6) / 2 / C DATA IMACH(7) / 2 / C DATA IMACH(8) / 15 / C DATA IMACH(9) / 32767 / C DATA IMACH(10)/ 2 / C DATA IMACH(11)/ 23 / C DATA IMACH(12)/ -128 / C DATA IMACH(13)/ 127 / C DATA IMACH(14)/ 39 / C DATA IMACH(15)/ -128 / C DATA IMACH(16)/ 127 / C C MACHINE CONSTANTS FOR THE HP 2100 C 4 WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA IMACH(1) / 5 / C DATA IMACH(2) / 6 / C DATA IMACH(3) / 4 / C DATA IMACH(4) / 1 / C DATA IMACH(5) / 16 / C DATA IMACH(6) / 2 / C DATA IMACH(7) / 2 / C DATA IMACH(8) / 15 / C DATA IMACH(9) / 32767 / C DATA IMACH(10)/ 2 / C DATA IMACH(11)/ 23 / C DATA IMACH(12)/ -128 / C DATA IMACH(13)/ 127 / C DATA IMACH(14)/ 55 / C DATA IMACH(15)/ -128 / C DATA IMACH(16)/ 127 / C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND C THE PERKIN ELMER (INTERDATA) 7/32. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 16 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / Z7FFFFFFF / C DATA IMACH(10) / 16 / C DATA IMACH(11) / 6 / C DATA IMACH(12) / -64 / C DATA IMACH(13) / 63 / C DATA IMACH(14) / 14 / C DATA IMACH(15) / -64 / C DATA IMACH(16) / 63 / C C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR). C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 5 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 5 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / "377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 54 / C DATA IMACH(15) / -101 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR). C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 5 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 5 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / "377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 62 / C DATA IMACH(15) / -128 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 32-BIT INTEGER ARITHMETIC. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 5 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 56 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 16-BIT INTEGER ARITHMETIC. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 5 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 16 / C DATA IMACH( 6) / 2 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 15 / C DATA IMACH( 9) / 32767 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 56 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILER C C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 1 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / O377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 60 / C DATA IMACH(15) /-1024 / C DATA IMACH(16) / 1023 / C C C MACHINE CONSTANTS FOR THE VAX 11/780 C DATA IMACH(1) / 5 / DATA IMACH(2) / 6 / DATA IMACH(3) / 5 / DATA IMACH(4) / 6 / DATA IMACH(5) / 32 / DATA IMACH(6) / 4 / DATA IMACH(7) / 2 / DATA IMACH(8) / 31 / DATA IMACH(9) /2147483647 / DATA IMACH(10)/ 2 / DATA IMACH(11)/ 24 / DATA IMACH(12)/ -127 / DATA IMACH(13)/ 127 / DATA IMACH(14)/ 54 / DATA IMACH(15)/ -127 / DATA IMACH(16)/ 127 / C C MACHINE CONSTANTS FOR THE Z80 MICROPROCESSOR C C DATA IMACH( 1) / 1/ C DATA IMACH( 2) / 1/ C DATA IMACH( 3) / 0/ C DATA IMACH( 4) / 1/ C DATA IMACH( 5) / 16/ C DATA IMACH( 6) / 2/ C DATA IMACH( 7) / 2/ C DATA IMACH( 8) / 15/ C DATA IMACH( 9) / 32767/ C DATA IMACH(10) / 2/ C DATA IMACH(11) / 24/ C DATA IMACH(12) / -127/ C DATA IMACH(13) / 127/ C DATA IMACH(14) / 56/ C DATA IMACH(15) / -127/ C DATA IMACH(16) / 127/ C C C***FIRST EXECUTABLE STATEMENT I1MACH IF (I .LT. 1 .OR. I .GT. 16) GO TO 10 C I1MACH=IMACH(I) RETURN C 10 CONTINUE WRITE(OUTPUT,9000) 9000 FORMAT('1ERROR 1 IN I1MACH - I OUT OF BOUNDS ') C C CALL FDUMP C C STOP END FUNCTION J4SAVE(IWHICH,IVALUE,ISET) C***BEGIN PROLOGUE J4SAVE C***DATE WRITTEN 790801 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. Z C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Saves and recalls some global variables needed by library C error handling routines C***DESCRIPTION C C Abstract C J4SAVE saves and recalls several global variables needed C by the library error handling routines. C C Description of Parameters C --Input-- C IWHICH - index of item desired. C = 1 refers to current error number. C = 2 refers to current error control flag. C = 3 refers to current unit number to which error C messages are to be sent. (0 means use standard.) C = 4 refers to the maximum number of times any C message is to be printed (as set by XERMAX). C = 5 refers to the total number of units to which C each error message is to be written. C = 6 refers to the 2nd unit for error messages C = 7 refers to the 3rd unit for error messages C = 8 refers to the 4th unit for error messages C = 9 refers to the 5th unit for error messages C IVALUE - the value to be set for the IWHICH-th parameter, C if ISET is .TRUE. . C ISET - if ISET=.TRUE., the IWHICH-th parameter will be C given the value, IVALUE. If ISET=.FALSE., the C IWHICH-th parameter will be unchanged, and IVALUE C is a dummy parameter. C --Output-- C The (old) value of the IWHICH-th parameter will be returned C in the function value, J4SAVE. C C Written by Ron Jones, with SLATEC Common Math Library Subcommittee C Adapted from Bell Laboratories PORT Library Error Handler C Latest revision --- 23 May 1979 C***REFERENCES JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR C HANDLING PACKAGE*, SAND78-1189, SANDIA LABORATORIES, C 1978. C***ROUTINES CALLED (NONE) C***END PROLOGUE J4SAVE C LOGICAL ISET INTEGER IPARAM(9) DATA IPARAM(1),IPARAM(2),IPARAM(3),IPARAM(4)/0,1,0,10/ DATA IPARAM(5)/1/ DATA IPARAM(6),IPARAM(7),IPARAM(8),IPARAM(9)/0,0,0,0/ C***FIRST EXECUTABLE STATEMENT J4SAVE J4SAVE = IPARAM(IWHICH) IF (ISET) IPARAM(IWHICH) = IVALUE RETURN END SUBROUTINE XERABT(MESSG,NMESSG) C***BEGIN PROLOGUE XERABT C***DATE WRITTEN 790801 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. R3C C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Aborts program execution and prints error message C***DESCRIPTION C C Abstract C ***Note*** machine dependent routine C XERABT aborts the execution of the program. C The error message causing the abort is given in the calling C sequence, in case one needs it for printing on a dayfile, C for example. C C Description of Parameters C MESSG and NMESSG are as in XERROR, except that NMESSG may C be zero, in which case no message is being supplied. C C Written by Ron Jones, with SLATEC Common Math Library Subcommittee C Latest revision --- 7 June 1978 C***REFERENCES JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR C HANDLING PACKAGE*, SAND78-1189, SANDIA LABORATORIES, C 1978. C***ROUTINES CALLED (NONE) C***END PROLOGUE XERABT C DIMENSION MESSG(NMESSG) C***FIRST EXECUTABLE STATEMENT XERABT STOP C RETURN END SUBROUTINE XERCTL(MESSG1,NMESSG,NERR,LEVEL,KONTRL) C***BEGIN PROLOGUE XERCTL C***DATE WRITTEN 790801 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. R3C C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Allows user control over handling of individual errors C***DESCRIPTION C C Abstract C Allows user control over handling of individual errors. C Just after each message is recorded, but before it is C processed any further (i.e., before it is printed or C a decision to abort is made), a call is made to XERCTL. C If the user has provided his own version of XERCTL, he C can then override the value of KONTROL used in processing C this message by redefining its value. C KONTRL may be set to any value from -2 to 2. C The meanings for KONTRL are the same as in XSETF, except C that the value of KONTRL changes only for this message. C If KONTRL is set to a value outside the range from -2 to 2, C it will be moved back into that range. C C Description of Parameters C C --Input-- C MESSG1 - the first word (only) of the error message. C NMESSG - same as in the call to XERROR or XERRWV. C NERR - same as in the call to XERROR or XERRWV. C LEVEL - same as in the call to XERROR or XERRWV. C KONTRL - the current value of the control flag as set C by a call to XSETF. C C --Output-- C KONTRL - the new value of KONTRL. If KONTRL is not C defined, it will remain at its original value. C This changed value of control affects only C the current occurrence of the current message. C***REFERENCES JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR C HANDLING PACKAGE*, SAND78-1189, SANDIA LABORATORIES, C 1978. C***ROUTINES CALLED (NONE) C***END PROLOGUE XERCTL C C***FIRST EXECUTABLE STATEMENT XERCTL RETURN END SUBROUTINE XERPRT(MESSG,NMESSG) C***BEGIN PROLOGUE XERPRT C***DATE WRITTEN 790801 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. Z C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Prints error messages C***DESCRIPTION C C Abstract C Print the Hollerith message in MESSG, of length NMESSG, C on each file indicated by XGETUA. C***REFERENCES JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR C HANDLING PACKAGE*, SAND78-1189, SANDIA LABORATORIES, C 1978. C***ROUTINES CALLED I1MACH,S88FMT,XGETUA C***END PROLOGUE XERPRT C INTEGER F(10),G(14),LUN(5) DIMENSION MESSG(NMESSG) DATA F(1),F(2),F(3),F(4),F(5),F(6),F(7),F(8),F(9),F(10) 1 / 1H( ,1H1 ,1HX ,1H, ,1H ,1H ,1HA ,1H ,1H ,1H) / DATA G(1),G(2),G(3),G(4),G(5),G(6),G(7),G(8),G(9),G(10) 1 / 1H( ,1H1 ,1HX ,1H ,1H ,1H ,1H ,1H ,1H ,1H / DATA G(11),G(12),G(13),G(14) 1 / 1H ,1H ,1H ,1H) / DATA LA/1HA/,LCOM/1H,/,LBLANK/1H / C PREPARE FORMAT FOR WHOLE LINES C***FIRST EXECUTABLE STATEMENT XERPRT NCHAR = I1MACH(6) NFIELD = 72/NCHAR CALL S88FMT(2,NFIELD,F(5)) CALL S88FMT(2,NCHAR,F(8)) C PREPARE FORMAT FOR LAST, PARTIAL LINE, IF NEEDED NCHARL = NFIELD*NCHAR NLINES = NMESSG/NCHARL NWORD = NLINES*NFIELD NCHREM = NMESSG - NLINES*NCHARL IF (NCHREM.LE.0) GO TO 40 DO 10 I=4,13 10 G(I) = LBLANK NFIELD = NCHREM/NCHAR IF (NFIELD.LE.0) GO TO 20 C PREPARE WHOLE WORD FIELDS G(4) = LCOM CALL S88FMT(2,NFIELD,G(5)) G(7) = LA CALL S88FMT(2,NCHAR,G(8)) 20 CONTINUE NCHLST = MOD(NCHREM,NCHAR) IF (NCHLST.LE.0) GO TO 30 C PREPARE PARTIAL WORD FIELD G(10) = LCOM G(11) = LA CALL S88FMT(2,NCHLST,G(12)) 30 CONTINUE 40 CONTINUE C PRINT THE MESSAGE NWORD1 = NWORD+1 NWORD2 = (NMESSG+NCHAR-1)/NCHAR CALL XGETUA(LUN,NUNIT) DO 50 KUNIT = 1,NUNIT IUNIT = LUN(KUNIT) IF (IUNIT.EQ.0) IUNIT = I1MACH(4) IF (NWORD.GT.0) WRITE (IUNIT,F) (MESSG(I),I=1,NWORD) IF (NCHREM.GT.0) WRITE (IUNIT,G) (MESSG(I),I=NWORD1,NWORD2) 50 CONTINUE RETURN END SUBROUTINE XERSAV(MESSG,NMESSG,NERR,LEVEL,ICOUNT) C***BEGIN PROLOGUE XERSAV C***DATE WRITTEN 800319 (YYMMDD) C***REVISION DATE 840405 (YYMMDD) C***CATEGORY NO. Z C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Records that an error occurred. C***DESCRIPTION C C Abstract C Record that this error occurred. C C Description of Parameters C --Input-- C MESSG, NMESSG, NERR, LEVEL are as in XERROR, C except that when NMESSG=0 the tables will be C dumped and cleared, and when NMESSG is less than zero the C tables will be dumped and not cleared. C --Output-- C ICOUNT will be the number of times this message has C been seen, or zero if the table has overflowed and C does not contain this message specifically. C When NMESSG=0, ICOUNT will not be altered. C C Written by Ron Jones, with SLATEC Common Math Library Subcommittee C Latest revision --- 19 Mar 1980 C***REFERENCES JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR C HANDLING PACKAGE*, SAND78-1189, SANDIA LABORATORIES, C 1978. C***ROUTINES CALLED I1MACH,S88FMT,XGETUA C***END PROLOGUE XERSAV C INTEGER F(17),LUN(5) DIMENSION MESSG(1) DIMENSION MESTAB(10),NERTAB(10),LEVTAB(10),KOUNT(10) DATA MESTAB(1),MESTAB(2),MESTAB(3),MESTAB(4),MESTAB(5), 1 MESTAB(6),MESTAB(7),MESTAB(8),MESTAB(9),MESTAB(10) 2 /0,0,0,0,0,0,0,0,0,0/ DATA NERTAB(1),NERTAB(2),NERTAB(3),NERTAB(4),NERTAB(5), 1 NERTAB(6),NERTAB(7),NERTAB(8),NERTAB(9),NERTAB(10) 2 /0,0,0,0,0,0,0,0,0,0/ DATA LEVTAB(1),LEVTAB(2),LEVTAB(3),LEVTAB(4),LEVTAB(5), 1 LEVTAB(6),LEVTAB(7),LEVTAB(8),LEVTAB(9),LEVTAB(10) 2 /0,0,0,0,0,0,0,0,0,0/ DATA KOUNT(1),KOUNT(2),KOUNT(3),KOUNT(4),KOUNT(5), 1 KOUNT(6),KOUNT(7),KOUNT(8),KOUNT(9),KOUNT(10) 2 /0,0,0,0,0,0,0,0,0,0/ DATA KOUNTX/0/ DATA F(1),F(2),F(3),F(4),F(5),F(6),F(7),F(8),F(9),F(10), 1 F(11),F(12),F(13),F(14),F(15),F(16),F(17) 2 /1H( ,1H1 ,1HX ,1H, ,1HA ,1H ,1H ,1H, ,1HI ,1H , 3 1H ,1H, ,1H2 ,1HI ,1H1 ,1H0 ,1H) / C***FIRST EXECUTABLE STATEMENT XERSAV IF (NMESSG.GT.0) GO TO 80 C DUMP THE TABLE IF (KOUNT(1).EQ.0) RETURN C PREPARE FORMAT NCHAR = I1MACH(6) CALL S88FMT(2,NCHAR,F(6)) NCOL = 20 - NCHAR CALL S88FMT(2,NCOL,F(10)) C PRINT TO EACH UNIT CALL XGETUA(LUN,NUNIT) DO 60 KUNIT=1,NUNIT IUNIT = LUN(KUNIT) IF (IUNIT.EQ.0) IUNIT = I1MACH(4) C PRINT TABLE HEADER WRITE (IUNIT,10) 10 FORMAT ('0 ERROR MESSAGE SUMMARY'/ 1 ' FIRST WORD NERR LEVEL COUNT') C PRINT BODY OF TABLE DO 20 I=1,10 IF (KOUNT(I).EQ.0) GO TO 30 WRITE (IUNIT,F) MESTAB(I),NERTAB(I),LEVTAB(I),KOUNT(I) 20 CONTINUE 30 CONTINUE C PRINT NUMBER OF OTHER ERRORS IF (KOUNTX.NE.0) WRITE (IUNIT,40) KOUNTX 40 FORMAT ('0OTHER ERRORS NOT INDIVIDUALLY TABULATED=',I10) WRITE (IUNIT,50) 50 FORMAT (1X) 60 CONTINUE IF (NMESSG.LT.0) RETURN C CLEAR THE ERROR TABLES DO 70 I=1,10 70 KOUNT(I) = 0 KOUNTX = 0 RETURN 80 CONTINUE C PROCESS A MESSAGE... C SEARCH FOR THIS MESSG, OR ELSE AN EMPTY SLOT FOR THIS MESSG, C OR ELSE DETERMINE THAT THE ERROR TABLE IS FULL. DO 90 I=1,10 II = I IF (KOUNT(I).EQ.0) GO TO 110 IF (MESSG(1).NE.MESTAB(I)) GO TO 90 IF (NERR.NE.NERTAB(I)) GO TO 90 IF (LEVEL.NE.LEVTAB(I)) GO TO 90 GO TO 100 90 CONTINUE C THREE POSSIBLE CASES... C TABLE IS FULL KOUNTX = KOUNTX+1 ICOUNT = 1 RETURN C MESSAGE FOUND IN TABLE 100 KOUNT(II) = KOUNT(II) + 1 ICOUNT = KOUNT(II) RETURN C EMPTY SLOT FOUND FOR NEW MESSAGE 110 MESTAB(II) = MESSG(1) NERTAB(II) = NERR LEVTAB(II) = LEVEL KOUNT(II) = 1 ICOUNT = 1 RETURN END SUBROUTINE XGETUA(IUNIT,N) C***BEGIN PROLOGUE XGETUA C***DATE WRITTEN 790801 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. R3C C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Returns unit number(s) to which error messages are being C sent C***DESCRIPTION C C Abstract C XGETUA may be called to determine the unit number or numbers C to which error messages are being sent. C These unit numbers may have been set by a call to XSETUN, C or a call to XSETUA, or may be a default value. C C Description of Parameters C --Output-- C IUNIT - an array of one to five unit numbers, depending C on the value of N. A value of zero refers to the C default unit, as defined by the I1MACH machine C constant routine. Only IUNIT(1),...,IUNIT(N) are C defined by XGETUA. The values of IUNIT(N+1),..., C IUNIT(5) are not defined (for N .LT. 5) or altered C in any way by XGETUA. C N - the number of units to which copies of the C error messages are being sent. N will be in the C range from 1 to 5. C C Written by Ron Jones, with SLATEC Common Math Library Subcommittee C***REFERENCES JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR C HANDLING PACKAGE*, SAND78-1189, SANDIA LABORATORIES, C 1978. C***ROUTINES CALLED J4SAVE C***END PROLOGUE XGETUA C DIMENSION IUNIT(5) C***FIRST EXECUTABLE STATEMENT XGETUA N = J4SAVE(5,0,.FALSE.) DO 30 I=1,N INDEX = I+4 IF (I.EQ.1) INDEX = 3 IUNIT(I) = J4SAVE(INDEX,0,.FALSE.) 30 CONTINUE RETURN END SUBROUTINE S88FMT(N,IVALUE,IFMT) C***BEGIN PROLOGUE S88FMT C***DATE WRITTEN 790801 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. Z C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Integer to character conversion C***DESCRIPTION C C Abstract C S88FMT replaces IFMT(1), ... ,IFMT(N) with the C characters corresponding to the N least significant C digits of IVALUE. C C Taken from the Bell Laboratories PORT Library Error Handler C Latest revision --- 7 June 1978 C***REFERENCES JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR C HANDLING PACKAGE*, SAND78-1189, SANDIA LABORATORIES, C 1978. C***ROUTINES CALLED (NONE) C***END PROLOGUE S88FMT C DIMENSION IFMT(N),IDIGIT(10) DATA IDIGIT(1),IDIGIT(2),IDIGIT(3),IDIGIT(4),IDIGIT(5), 1 IDIGIT(6),IDIGIT(7),IDIGIT(8),IDIGIT(9),IDIGIT(10) 2 /1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9/ C***FIRST EXECUTABLE STATEMENT S88FMT NT = N IT = IVALUE 10 IF (NT .EQ. 0) RETURN INDEX = MOD(IT,10) IFMT(NT) = IDIGIT(INDEX+1) IT = IT/10 NT = NT - 1 GO TO 10 END