PROGRAM MAINPE
C
C     ------------------------------------------------------------
C
C     Testexample for PERIOD - Lorenz-Rayleigh-Benard-Problem
C
C*  Purpose           Testexample for code PERIOD
C*  Version           0.5 (Test Version)
C*  Revision          July 1984
C*  Latest Change     January 1991
C*  Library           CodeLib
C*  Code              Fortran 77, Double Precision
C*  Environment       Standard Fortran 77 environment on PC's,
C                     workstations and hosts.
C
C     ------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION Z(3),Y(3,11),X(11),RW(6000),IW(100),ICALL(8),USCAL(3)
     1          ,FM(3,2),IPAR(10),RH(2000),IH(100)
      REAL CPS,XZEIT
      INTEGER IFAIL
      EXTERNAL PRDIFX,FCN
      COMMON /DATA/SIGMA,R,B
C
C  NUMBER OF ORDINARY DIFFERENTIAL EQUATIONS
      N=3
C  NUMBER OF SHOOTING NODES
      M=11
C
      M1=M-1
C
C  AMOUNT OF REAL- AND INTEGER-WORKSPACE
      NRW=6000
      NIW=100
C
C  SHOOTING NODES
      X(1)=0.D0
      RED=1.D0/DBLE(M1)
      DO 1 I=2,M
1     X(I)=X(I-1)+RED
C
      SIGMA=1.6D1
      R=153.083D0
      B=4.D0
C
C  ESTIMATE OF PERIOD
      P=0.95D0
C
C  INITIAL VALUES OF TRAJECTORY
      Y(1,1)=0.D0
      Y(2,1)=-28.D0
      Y(3,1)=140.D0
C
C  CALLING DIFEXP TO GENERATE INITIAL TRAJECTORY
      NRH=2000
      NIH=100
      IPAR(1)=0
      IPAR(2)=0
      IPAR(3)=0
      IPAR(4)=N
      IPAR(5)=0
      IPAR(6)=0
      T=0.D0
      H=1.D-3
      Z(1)=Y(1,1)
      Z(2)=Y(2,1)
      Z(3)=Y(3,1)
      TOL=1.D-3
      DO 22 I=2,M
      TEND=X(I)*P
      HMAX=TEND-T
      CALL PRDIFX (N,FCN,T,Z,TEND,TOL,HMAX,H,HS,USCAL,NRH,RH,NIH,IH,
     2             IPAR)
      Y(1,I)=Z(1)
      Y(2,I)=Z(2)
      Y(3,I)=Z(3)
22    CONTINUE
C
C
C  DESIRED RELATIVE ACCURACY FOR SOLUTION
      EPS=1.D-6
C  CLASSIFICATION OF RIGHT-HAND SIDE (AUTONOMOUS SYSTEM)
      ICALL(1)=1
C  PROBLEM IS HIGHLY NONLINEAR
      ICALL(2)=2
C  MAXIMUM PERMITTED NEWTON ITERATIONS
      ICALL(3)=30
C  SOLVING VARIATIONAL EQUATION
      ICALL(4)=1
C  RANK-1 UPDATES ALLOWED
      ICALL(5)=1
C  ITERATIVE REFINEMENT IS ACTIVATED
      ICALL(6)=1
C  THE SYSTEM IS NONSTIFF
      ICALL(7)=0
C  PRINT PARAMETER
      ICALL(8)=0
C
      CALL ZIBSEC(CPS,IFAIL)
C
      CALL PERIOD(PRDIFX,N,M,X,Y,P,EPS,FM,6000,RW,100,IW,ICALL)
C
      CALL ZIBSEC(XZEIT,IFAIL)
      XZEIT=XZEIT-CPS
      WRITE(6,666) XZEIT
666   FORMAT(1X,/,' CPU-TIME=',F5.2)
      STOP
      END
C
      SUBROUTINE FCN(T,Z,DZ)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Z(3),DZ(3)
      COMMON /DATA/SIGMA,R,B
      DZ(1)=SIGMA*(Z(2)-Z(1))
      DZ(2)=Z(1)*(R-Z(3))-Z(2)
      DZ(3)=Z(1)*Z(2)-B*Z(3)
      RETURN
      END
C
      SUBROUTINE DFDY(T,Z,DF)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Z(3),DF(3,3)
      COMMON /DATA/SIGMA,R,B
      DATA ZERO/0.D0/,ONE/1.D0/
      DF(1,1)=-SIGMA
      DF(1,2)=SIGMA
      DF(1,3)=ZERO
      DF(2,1)=R-Z(3)
      DF(2,2)=-ONE
      DF(2,3)=-Z(1)
      DF(3,1)=Z(2)
      DF(3,2)=Z(1)
      DF(3,3)=-B
      RETURN
      END
      REAL FUNCTION SECOND(IDUMMY)
C
C     This subroutine is machine dependent. It's contents
C     may be replaced by code, which delivers the actual
C     time stamp, to get a time measurement. 
C
      SECOND=0.0E0
      RETURN
      END