PROGRAM PLFADR IMPLICIT DOUBLE PRECISION (A-H, O-Z) C C C* Begin Prologue PLFADR C C --------------------------------------------------------------------- C C* Purpose Driver testing MEXX with C example pendulum (PL) C using full (F) matrix mode C with options (A) in order to get an integration C monitor and solution output after each integration C step C C* File plfadr.f C* Version 1.1 C* Latest Change 96/07/19 (1.6) C* Modification history C 1.1 Correct workspace formula. C Avoid free format because of compiler dependency. C* Library CodeLib C* Code Fortran 77 C Double Precision C* Environment Standard version for FORTRAN77 environments on C PCs, workstations, and hosts C C* Subroutines supplied C C PLFIN Provide initial values C PLFF Specify the problem C PLOUT Routine for solution output C DUDOUT Dummy routine (instaed of dense output) C DUSWI Dummy routine (instead of switching functions) C C Problem dimensions for a maximum of 30 bodies. C The actual number of bodies will be requested from standard input. C PARAMETER ( $ NBODY = 30, $ NPDIM = 3*NBODY, $ NVDIM = NPDIM, $ NUDIM = 0, $ NGDIM = 2*NBODY, $ NLDIM = NGDIM, $ NPVUD = NPDIM*2 + NUDIM, $ NZGDIM = 8*NBODY - 4 $ ) C INTEGER NL1, NU1, NG1, NGL C PARAMETER ( NL1 = MAX (1,NLDIM) , C $ NU1 = MAX (1,NUDIM) , C $ NG1 = MAX (1,NGDIM) , C $ NGL = MAX (1,NGDIM,NLDIM) ) C C Workspace dimension C PARAMETER ( $ LIWK = NPDIM + 4*NVDIM + NLDIM + NUDIM + 60, $ LRWK = (NVDIM+NLDIM+1)**2 + NPDIM*(NGDIM+NLDIM+1+18) + $ NVDIM*(NVDIM+45) + 28*(NLDIM+1) + 18*(NUDIM+1) + $ NGDIM+1 + 50 $ ) C PARAMETER ( C $ LIWK = NPDIM + 4*NVDIM + NLDIM + NUDIM + 60, C $ LRWK = (NVDIM+NL1)**2 + NPDIM*(NGL+18) + C $ NVDIM*(NVDIM+45) + 28*NL1 + 18*NU1 + NG1 + 50 C $ ) C EXTERNAL PLFIN, PLFF, PLOUT, DUDOUT, DUSWI C INTEGER MXJOB(150), IWK(LIWK) DOUBLE PRECISION RWK(LRWK) DIMENSION P(NPDIM), V(NVDIM), RTOL(NPVUD), ATOL(NPVUD) DIMENSION A(NVDIM), RLAM(NLDIM), U(NUDIM+1) C COMMON / PEND / PMASS(NBODY), PLEN(NBODY), NPEND C CHARACTER*30 PROTXT DATA PROTXT /'Pendulum'/ C CHARACTER*67 CMODUL DATA CMODUL /'@(#) plfadr.f 1.6 from 96/07/19 (mexx 1.2) $ '/ C DATA MXJOB /150*0/, RWK /LRWK*0.0D0/, IWK /LIWK*0/ C C -- Machine constants LUSIN = 5 LUSOUT = 6 C C -- Start message WRITE (LUSOUT,9000) ' ' WRITE (LUSOUT,9000) $ ' ====================================================' WRITE (LUSOUT,9000) CMODUL WRITE (LUSOUT,9000) $ ' ====================================================' WRITE (LUSOUT,9000) ' ' C C -- Initialize problem NP = NPDIM NG = NGDIM NU = NUDIM NV=NP NL=0 NBLK=0 NMRC=0 NZGMAX=0 NZFMAX=0 NSWMAX=0 C CALL PLFIN (NP, NV, NL, NG, NU, $ T0, P, V, U, A, RLAM, PROTXT) C C Set options for integration C T=T0 TEND = 1.D0 H = 1.D-3 ITOL = 0 TOL = 1.D-3 RTOL(1) = TOL ATOL(1) = TOL C C Select I/O options C - General print output MXJOB(11) = 1 MXJOB(12) = LUSOUT C - Integration monitor C Note step size reductions MXJOB(13) = 2 MXJOB(14) = LUSOUT C C Select options for solution output C MOUT = 1 C MXJOB(30) = MOUT C IDID = 0 C CALL MEXX (NP, NV, NL, NG, NU, $ PLFF, T, TEND, P, V, U, A, RLAM, $ ITOL, RTOL, ATOL, $ H, MXJOB, IDID, LIWK, IWK, LRWK, RWK, $ PLOUT, DUDOUT, DUSWI) C IF (IDID .EQ. 0) THEN WRITE (LUSOUT,9000) ' o.k.' ELSE WRITE (LUSOUT,9100) ' Error return code from MEXX:', IDID ENDIF C C WRITE (LUSOUT,9000) ' ' WRITE (LUSOUT,9200) ' Position values at the end point:', T WRITE (LUSOUT,9000) ' ' WRITE (LUSOUT,9000) ' P(I),I=1,2,..' WRITE (LUSOUT,9000) ' ' WRITE (LUSOUT,9210) (P(I),I=1,NP) C C WRITE (LUSOUT,9100) ' ' WRITE (LUSOUT,9100) 'C===================================' WRITE (LUSOUT,9100) ' ' WRITE (LUSOUT,9200) ' TOL ', RTOL(1) WRITE (LUSOUT,9100) ' ' WRITE (LUSOUT,9100) ' Statistics: NSTEP, NACCPT, NREJCT' WRITE (LUSOUT,9100) ' ', MXJOB(51), MXJOB(52), MXJOB(53) WRITE (LUSOUT,9100) ' ' WRITE (LUSOUT,9100) ' NFPROB, NFCN, NGCN' WRITE (LUSOUT,9100) ' ', MXJOB(54), MXJOB(55), MXJOB(56) WRITE (LUSOUT,9100) ' ' WRITE (LUSOUT,9100) ' NMULT, NDEC, NSOL' WRITE (LUSOUT,9100) ' ', MXJOB(57), MXJOB(58), MXJOB(59) WRITE (LUSOUT,9100) ' ' WRITE (LUSOUT,9100) 'C===================================' C C WRITE (LUSOUT,9000) ' ' WRITE (LUSOUT,9000) ' End of example ', PROTXT WRITE (LUSOUT,9000) ' Bye' C STOP C 9000 FORMAT (2A) 9100 FORMAT (A, 3I8) 9200 FORMAT (A, 1PD8.1) 9210 FORMAT ((4(1PD19.10))) C C end of program PLFADR C END C SUBROUTINE PLFIN (NP, NV, NL, NG, NU, $ T0, P, V, U, A, RLAM, PROTXT) C IMPLICIT DOUBLE PRECISION (A-H, O-Z) C INTEGER NP, NV, NL, NG, NU DIMENSION P(*), V(*), U(*), A(*), RLAM(*) CHARACTER*(*) PROTXT C C ********************************************************************** C C Pendulum (variable number of bodies < 31) C C C Parameter: masses pmass(i) C lengths plen(i) C C C ********************************************************************** C C* File plfin.f C* Version 1.1 C* Latest Change 96/07/19 (1.4) C* Modification history C 1.1 Correct initial values of RLAM. C New free format starts text output with ' '. C C ------------------------------------------------------- C PARAMETER (NBMAX = 30) COMMON / PEND / PMASS(NBMAX), PLEN(NBMAX), NPEND C SAVE C CHARACTER*67 CMODUL DATA CMODUL /'@(#) plfin.f 1.4 from 96/07/19 (mexx 1.2) $ '/ C PROTXT = 'Pendulum (default options)' C C Machine constants (logical units stdin and stdout) LUSIN = 5 LUSOUT = 6 C WRITE (LUSOUT,*) 'Start example ', PROTXT WRITE (LUSOUT,*) ' ' C C Dimensions C Get problem parameters C WRITE (LUSOUT,*) 'Enter number of bodies (default 8):' C READ (LUSIN,*, ERR=1000, END=1000) NPEND C GOTO 1010 1000 NPEND = 8 1010 WRITE (LUSOUT,*) 'Given:', NPEND IF (NPEND .LT. 1 .OR. NPEND .GT. NBMAX) THEN WRITE (LUSOUT,*) $ 'Number of bodies is out of range (default value chosen)' NPEND = 8 ENDIF WRITE (LUSOUT,*) ' ' C NP = 3*NPEND NV = NP NG = 2*NPEND NL = NG NU = 0 C C Model parameters T0 = 0.D0 PMASSH = 1.0D0/DBLE(NPEND) PLENH = 0.5D0/DBLE(NPEND) DO 1020 I=1,NPEND PMASS(I) = PMASSH PLEN(I) = PLENH 1020 CONTINUE C C Initial values PIHALF = ASIN(1.D0) P(1) = PLEN(1) P(2) = 0.D0 P(3) = PIHALF C V(1) = 0.D0 V(2) = 0.D0 V(3) = 0.D0 C RLAM(1) = 0.D0 RLAM(2) = 0.D0 C DO 1030 I=2,NPEND L = (I-1)*3 + 1 P(L) = P(L-3) + PLEN(I-1) + PLEN(I) P(L+1) = P(2) P(L+2) = P(3) V(L) = 0.D0 V(L+1) = 0.D0 V(L+2) = 0.D0 L = (I-1)*2 + 1 RLAM(L) = 0.D0 RLAM(L+1) = 0.D0 1030 CONTINUE C RETURN END C SUBROUTINE PLFF (NP, NV, NL, NG, NU, LDG, $ T, P, V, U, RLAM, $ AM, GP, F, PDOT, UDOT, G, GI, FL, $ QFLAG, IFAIL) C IMPLICIT DOUBLE PRECISION (A-H, O-Z) C Input parameters: INTEGER NP, NV, NL, NG, NU, LDG LOGICAL QFLAG(9) DOUBLE PRECISION T, P(NP), V(NV), U(NU), RLAM(NL) C Output parameters: DOUBLE PRECISION AM(NV,NV), GP(LDG,*), F(NV), PDOT(NP), UDOT(NU), $ G(NG), GI(NG), FL(NV,NL) C C C* File plff.f C* Version 1.1 C* Latest Change 96/05/03 (1.5) C* Modification history C 1.1 New free format starts text output with ' '. C C PARAMETER (NBMAX = 30) COMMON / PEND / PMASS(NBMAX), PLEN(NBMAX), NPEND SAVE C DATA GRAV/10./ C C Machine constant (logical unit stdout) LUSOUT = 6 C IFAIL = 0 C IF (QFLAG(1)) THEN DO 1010 J=1,NV DO 1000 I=1,NV AM(I,J) = 0.D0 1000 CONTINUE 1010 CONTINUE DO 1020 I=1,NPEND L = 3*(I-1) + 1 AM(L ,L ) = PMASS(I) AM(L+1,L+1) = PMASS(I) AM(L+2,L+2) = PLEN(I)**2*PMASS(I)/3.D0 1020 CONTINUE ENDIF C IF (QFLAG(2) .OR. QFLAG(3)) THEN DO 1040 J=1,NV DO 1030 I=1,NG GP(I,J) = 0.D0 1030 CONTINUE 1040 CONTINUE DO 1050 I=1,NPEND IROW = (I-1)*2 + 1 ICOL = (I-1)*3 + 1 GP(IROW, ICOL) = -1.D0 GP(IROW, ICOL+2) = PLEN(I)*COS(P(ICOL+2)) GP(IROW+1, ICOL+1) = -1.D0 GP(IROW+1, ICOL+2) = PLEN(I)*SIN(P(ICOL+2)) IF (I .GT. 1) THEN GP(IROW, ICOL-3) = 1.D0 GP(IROW, ICOL-1) = PLEN(I-1)*COS(P(ICOL-1)) GP(IROW+1, ICOL-2) = 1.D0 GP(IROW+1, ICOL-1) = PLEN(I-1)*SIN(P(ICOL-1)) ENDIF 1050 CONTINUE ENDIF IF (QFLAG(4)) THEN DO 1060 I=1,NPEND L = 3*(I-1) + 1 F(L)= 0.D0 F(L+1) = -GRAV*PMASS(I) F(L+2) = 0.D0 1060 CONTINUE ENDIF IF (QFLAG(5)) THEN DO 1070 I=1,NV PDOT(I) = V(I) 1070 CONTINUE ENDIF IF (QFLAG(7)) THEN G(1) = -P(1) + PLEN(1)*SIN(P(3)) G(2) = -P(2) - PLEN(1)*COS(P(3)) DO 1080 I=2,NPEND L2 = 2*(I-1) + 1 L3 = 3*(I-1) + 1 G(L2) = P(L3-3) + PLEN(I-1)*SIN(P(L3-1)) $ -P(L3) + PLEN(I)*SIN(P(L3+2)) G(L2+1) = P(L3-2) - PLEN(I-1)*COS(P(L3-1)) $ -P(L3+1) - PLEN(I)*COS(P(L3+2)) 1080 CONTINUE ENDIF IF (QFLAG(8)) THEN DO 1090 I=1,NG GI(I) = 0.D0 1090 CONTINUE ENDIF IF (QFLAG(9)) THEN DO 1110 J=1,NL DO 1100 I=1,NV FL(I,J) = 0.D0 1100 CONTINUE 1110 CONTINUE WRITE (LUSOUT,*) 'PLFF (FL not available)' ENDIF C RETURN END C SUBROUTINE PLOUT (NP, NV, NU, NL, T, P, V, U, A, RL, INFO, IRTRN) C IMPLICIT DOUBLE PRECISION (A-H, O-Z) C INTEGER INFO(50), NP, NV, NU INTEGER IRTRN DOUBLE PRECISION T DIMENSION P(NP), V(NV), U(NU), A(NV), RL(NL) C C* Begin Prologue PLOUT C C* Purpose Write position values at integration points of MEXX C C* File plout.f C* Version 1.0 C* Latest Change 93/09/15 (1.1) C C C* Input arguments: C C NP Int Length of the position vector P C C NV Int Length of the velocity vector V C C NU Int Length of the vector of dynamic variables U C C NA Int Length of the acceleration vector A C C NL Int Length of the vector of Lagrange multipliers C C T Dble Current time C C P(NP) Dble The position vector. C C V(NV) Dble The velocity vector. C C U(NU) Dble The dynamic variable vector. C C A(NV Dble the dense output acceleration vector. C C RL(NL) Dble the dense output Lagrange multiplier vector. C C INFO(50) Int Info array C INFO(1) .. INFO(40) see MXJOB(1) .. MXJOB(40) C in MEXX C INFO(41): Type of call: C = 0: the first call C = 1: an intermediate call C = 2: the last call at TFIN C = 3: at a root of switching function C and integration continues C = 4: at a root of switching function C and integration will terminate C (this is the last call) C C INFO(42): number of current integration step C C INFO(43): interpolation indicator C = 0: current t is an integration point C = 1: the solution values are interpolated C values C C INFO(44): call counter C C Do NOT alter any of the input arguments C C* Output argument: C C IRTRN Int Stop indicator C If IRTRN .NE. 0 on return, MEXX will terminate C C SAVE C C Logical units C LUSOL output of solutions C DATA LUSOL /14/ C IRTRN = 0 IF (INFO(41) .EQ. 0) THEN OPEN (LUSOL, FILE = 'PLOUT.SOL') ENDIF C C WRITE (LUSOL,9000) INFO(42), T, (P(I),I=1,NP) C IRTRN = 0 C RETURN 9000 FORMAT (I19, 3(1PD19.10), /, (4(1PD19.10))) END C SUBROUTINE DUDOUT (NDP, NDV, NDU, NDA, NDL, T, $ DP, DV, DU, DA, DL, INDP, $ INDV, INDU, INDA, INDL, INFO, IFAIL) C IMPLICIT DOUBLE PRECISION (A-H, O-Z) C INTEGER NDP, NDV, NDU, NDA, NDL DOUBLE PRECISION T, DP(*), DV(*), DU(*), DA(*), DL(*) INTEGER INDP(*), INDV(*), INDU(*), INDA(*), INDL(*), INFO(50), $ IFAIL C C C* File dudout.f C* Version 1.0 C* Latest Change 96/03/05 (1.2) C C* Input arguments: C C NDP Int Number of dense output components of position C vector P. Dimension of arrays DP and INDP. The C current value is either the user prescibed value C IWK(51) (if MDIND=INFO(33)=1 was set) or NP, the C dimension of the position vector P (set by MEXX if C INFO(33)=0 was set). C C NDV Int as NDP but now for the velocity vector V. C C NDU Int as NDP but now for the vector of dynamic C variables U. C C NDA Int as NDP but now for the acceleration vector A. C C NDL Int as NDP but now for the vector of Lagrange C multipliers. C C T Dble Current time C C DP Dble the dense output position vector. C C DV Dble the dense output velocity vector. C C DU Dble the dense output dynamic variable vector. C C DA Dble the dense output acceleration vector. C C DL Dble the dense output Lagrange multiplier vector. C C INDP(NDP) Int indicator array of length NP. List of components of C the position vecor p which was selected for dense C output. INDP(I)=J means that component J of the C vector P is available as component I of vector DP. C (I.e. DP(I)=P(INDP(I)) C Note, that MEXX automatically provides INDP(I)=I C if INFO(33)=0 was set. C C INDV(NDV) Int as INDP but now for the velocity vector V. C C INDU(NDU) Int as INDP but now the vector of dynamic variables U. C C INDU(NDA) Int as INDP but now for the acceleration vector A. C C INDL(NDL) Int as INDP but now the vector of Lagrange multipliers. C C INFO(50) Int Info array C INFO(1) .. INFO(40) see MXJOB(1) .. MXJOB(40) C in MEXX C INFO(41): Type of call: C = 0: the first call C = 1: an intermediate call C = 2: the last call at TFIN C = 3: at a root of switching function C and integration continues C = 4: at a root of switching function C and integration will terminate C (this is the last call) C C INFO(42): number of current integration step C C INFO(43): interpolation indicator C = 0: current t is an integration point C = 1: the solution values are interpolated C values C C INFO(44): call counter C C IFAIL Int Error indicator. IFAIL=0 on input C C Do NOT alter any of the input arguments. C C C* Output arguments: C C IFAIL Int Error indicator C If IFAIL.LT.0 on return, MEXX will terminate C C---------------------------------------------------------------------- C C LUSOUT = 6 C WRITE (LUSOUT,*) ' +++ Warning in DUDOUT' WRITE (LUSOUT,*) $ ' This dummy routine should never have been called.' C IFAIL = 0 C RETURN C C end of subroutine DUDOUT C END C SUBROUTINE DUSWI (NP, NV, NL, NU, T, P, V, U, A, RLAM, NSWIF, G) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C DIMENSION P(NP), V(NV), U(*), A(NV), RLAM(NL), G(*) C C* File duswi.f C* Version 1.0 C* Latest Change 96/03/05 (1.2) C C LUSOUT = 6 C WRITE (LUSOUT,*) ' +++ Warning in DUSWI' WRITE (LUSOUT,*) $ ' This dummy routine should never have been called.' C IFAIL = 0 C RETURN C END