C * * * * * * * * * * * * * * * * * * * * * * * * *
C --- DRIVER FOR RADAU5: FRACTIONAL PROBLEM IN DIETHELM ET AL.
C     version of 27 May 2025
C * * * * * * * * * * * * * * * * * * * * * * * * *
        include '../../radau5.f'
        include '../../decsol.f'
        include 'dc_sumexp.f'
        IMPLICIT REAL*8 (A-H,O-Z)
C --- PARAMETERS FOR RADAU5
        INTEGER, PARAMETER :: NDM=5000
        PARAMETER (LWORK=4*NDM*NDM+12*NDM+20,LIWORK=3*NDM+20)               
        REAL*8, dimension(20) :: RPAR
        INTEGER, dimension(20) :: IPAR
        DIMENSION Y(NDM),WORK(LWORK),IWORK(LIWORK),YEX(NDM)
        DIMENSION ATOL(NDM),RTOL(NDM)
        EXTERNAL FGAMI,JGAMI,JGAMIEX,MAS,SOLOUT
        REAL*4 start,finish
C --- PARAMETERS OF THE PROBLEM
        alpha=0.5D0
        RPAR(1)=alpha
C --- DESIRED ACCURACY
        TOL=1.0D-8
        EPSIL=TOL
C --- END OF INTEGRATION
        XEND=1.0D0
        CALL HMNEXP (XEND,ALPHA,EPSIL,DELTA,HEXP,MEXP,NEXP)
        write (6,*) '  '
        WRITE(*,*) 'DELTA,HEXP,MEXP,NEXP ', delta, HEXP, MEXP, NEXP
        RPAR(6)=HEXP
        IPAR(1)=MEXP
        IPAR(2)=NEXP
C --- SAME TOL FOR ALL COMNPONENTS (ITOL = 0) COMPONENTWISE TOL (ITOL = 1)
        DO ITOL = 0,0
C --- FULL LINEAR ALGEBRA (ILIN = 1) AND ADAPTED TO SUMEXP (ILIN = 2)
        DO ILIN = 1,2
           write (6,*) '  '
           MI=1
           ND=1
C --- ND=1 equation and 1 integral (as y-component)
           N=ND+1+MI*(NEXP-MEXP)
           IF (ILIN.EQ.2) N=N+2
C --- INITIAL VALUES
           X=0.0D0 
           Y=0.0D0
           Y(1)=0.D0
C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM? NO
           IMAS=1
           MLMAS=0
           MUMAS=0
C --- OUTPUT ROUTINE IS USED DURING INTEGRATION? NO
           IOUT=0
c           IOUT=1
C --- REQUIRED TOLERANCE
           RTOL=TOL
           ATOL=TOL
           IF (ITOL.EQ.1) THEN
             FACTOL = SQRT(TOL)/10.0D0
             DO I=ND+1,N
               ATOL(I)=FACTOL
               RTOL(I)=FACTOL
             END DO
           END IF
           H=0.001D0
C --- SET DEFAULT VALUES 
           DO I=1,20
             IWORK(I)=0
             WORK(I)=0.D0
           END DO
C ==================================================================
C          FULL LIN. ALGEBRA OPTION           
C ==================================================================           
           IF (ILIN.EQ.1) THEN
             IJAC=1
c             IJAC=0
             MLJAC=N
C --- CALL OF THE SUBROUTINE RADAU5
             CALL cpu_time(start)
             CALL RADAU5(N,FGAMI,X,Y,XEND,H,
     &                  RTOL,ATOL,ITOL,
     &                  JGAMI,IJAC,MLJAC,MUJAC,
     &                  MAS,IMAS,MLMAS,MUMAS,
     &                  SOLOUT,IOUT,
     &                  WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID)
             CALL cpu_time(finish)
           END IF
C ==================================================================
C          SUMEXP LIN. ALGEBRA OPTION           
C ==================================================================           
           IF (ILIN.EQ.2) THEN
             IJAC=1
             MUJAC=0
             MLJAC=2
C --- CALL OF THE SUBROUTINE RADAU5
             CALL cpu_time(start)
             CALL RADAU5(N,FGAMI,X,Y,XEND,H,
     &                  RTOL,ATOL,ITOL,
     &                  JGAMIEX,IJAC,MLJAC,MUJAC,
     &                  MAS,IMAS,MLMAS,MUMAS,
     &                  SOLOUT,IOUT,
     &                  WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID)
             CALL cpu_time(finish)
           END IF
           YEX(1)=(XEND**4-1.5D0*XEND**(alpha/2))**2
           ERR=0.0D0
           DO I=1,1
             ERR = ERR + (Y(I)-YEX(I))**2
           END DO
           ERR=SQRT(ERR)/YEX(1)
C --- PRINT FINAL SOLUTION
           write (6,77) ILIN,TOL,finish-start
 77        FORMAT(1X,'radau5, lin-alg = ',I1,'   TOL = ',E9.3,
     &           '    time = ',E8.2)
           WRITE (6,99) X,Y(1),ERR
 99        FORMAT(1X,'X =',F10.3,'    Y =',E20.12,'    ERR = ',E9.3)
           WRITE (6,91) (IWORK(J),J=14,20)
 91        FORMAT(' fcn=',I5,' jac=',I4,' step=',I4,' accpt=',I4,
     &        ' rejct=',I3,' dec=',I4,' sol=',I5)
        END DO
        END DO
        STOP
        END
C
        SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN)
C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS BY USING "CONTR5"
        IMPLICIT REAL*8 (A-H,O-Z)
        REAL*8, dimension(20) :: RPAR
        INTEGER, dimension(20) :: IPAR
        DIMENSION Y(N),CONT(LRC)
        COMMON /INTERN/XOUT
        IF (NR.EQ.1) THEN
           XOUT=0.0D0
        ELSE
 10        CONTINUE
           IF (X.GE.XOUT) THEN
              WRITE (6,99) XOUT,CONTR5(1,XOUT,CONT,LRC), NR-1        
              XOUT=XOUT+1.0D-2
              GOTO 10
           END IF
        END IF
 99     FORMAT(1X,'X =',F14.10,'    Y =',E18.10,'    NSTEP = ',I4)
        RETURN
        END
C
        SUBROUTINE FGAMI(N,X,Y,F,RPAR,IPAR)
C --- RIGHT-HAND SIDE OF EQUATION
        IMPLICIT REAL*8 (A-H,O-Z)
        REAL*8, dimension(20) :: RPAR
        INTEGER, dimension(20) :: IPAR
        DIMENSION Y(N),F(N)
C       Initial condition  
        Y(1)=MAX(0.0D0,Y(1))      
        Y0 = 0.D0
        alpha = RPAR(1)
        MEXP=IPAR(1)
        NEXP=IPAR(2)
        LEXP=NEXP-MEXP
        HEXP=RPAR(6)    
        PI=4.0D0*ATAN(1.0D0)
        F1 = (Gamma(9.0D0)/Gamma(9.D0-alpha))*X**(8.D0-alpha) - 3*
     1     Gamma(5.D0+alpha/2)/Gamma(5.D0-alpha/2)*X**(4.D0-alpha/2)+
     2     9.D0/4.D0*Gamma(alpha+1.D0)+(3.D0/2.D0*X**(alpha/2)-X**4)**3
     3     -Y(1)**(3.D0/2.D0) 
        APPINT=0.0D0
        DO I=3,LEXP+2
           IND=MEXP+(I-3)
           GAMMAI =  EXP(IND*HEXP) 
           CI     =  EXP((1-alpha)*IND*HEXP)
           F(I)   = - GAMMAI*Y(I) + F1
           APPINT =  APPINT + CI*Y(I)
        END DO
        APPINT=APPINT*HEXP*SIN(ALPHA*PI)/PI
        F(1) = Y0 + Y(2) - Y(1)
        F(2) = APPINT - Y(2)
        RETURN
        END
C
        SUBROUTINE JGAMI(N,X,Y,DFY,LDFY,RPAR,IPAR)
C --- JACOBIAN with full linear algebra
        IMPLICIT REAL*8 (A-H,O-Z)
        REAL*8, dimension(20) :: RPAR
        INTEGER, dimension(20) :: IPAR
        DIMENSION Y(N),DFY(LDFY,N)
        Y(1)=MAX(0.0D0,Y(1))      
        alpha = RPAR(1)
        MEXP=IPAR(1)
        NEXP=IPAR(2)
        LEXP=NEXP-MEXP
        HEXP=RPAR(6)    
        DFY = 0.0D0    
        DFY(1,1) = -1.0D0
        DFY(1,2) =  1.0D0
        DFY(2,2) = -1.0D0
        PI=4.0D0*ATAN(1.0D0)
        FAC=HEXP*SIN(ALPHA*PI)/PI
        II=3
        DO I=1,LEXP
          IND=MEXP+(I-1)
          GAMMAI =  EXP(IND*HEXP) 
          CI =      EXP((1-alpha)*IND*HEXP)
          DFY(2,II) = CI*FAC
          DFY(II,II) = - GAMMAI
          DFY(II,1) = -(3.D0/2.D0)*Y(1)**(1.D0/2.D0)
          II=II+1
        END DO
        RETURN
        END
C
C ==================================================================
C          SUMEXP LIN. ALGEBRA OPTION           
C ==================================================================    
        SUBROUTINE JGAMIEX(N,X,Y,DFY,LDFY,RPAR,IPAR)
C --- JACOBIAN for the option "sumexp" in the linear algebra
        IMPLICIT REAL*8 (A-H,O-Z)
        REAL*8, dimension(20) :: RPAR
        INTEGER, dimension(20) :: IPAR       
        DIMENSION Y(N),DFY(LDFY,N)
        Y(1)=MAX(0.0D0,Y(1))      
        alpha = RPAR(1)
        MEXP=IPAR(1)
        NEXP=IPAR(2)
        LEXP=NEXP-MEXP
        HEXP=RPAR(6)        
        DFY = 0.0D0    
        DFY(1,1)=-1.D0
        DFY(1,2)=1.0D0
        DFY(2,2)=-1.0D0
        DFY(2,3)=1.0D0
        DFY(1,4)=-1.5D0*Y(1)**(0.5D0)
        PI=4.0D0*ATAN(1.0D0)
        FAC=HEXP*SIN(ALPHA*PI)/PI
        II=5
        DO I=1,LEXP
          IND=MEXP+(I-1)
          GAMMAI =  EXP(IND*HEXP) 
          CI =      EXP((1-alpha)*IND*HEXP)
          DFY(1,II) = CI*FAC
          DFY(2,II) = - GAMMAI
          II=II+1
        END DO
        RETURN
        END
C
         SUBROUTINE MAS(N,AM,LMAS,RPAR,IPAR)
         IMPLICIT REAL*8 (A-H,O-Z)
         DIMENSION AM(LMAS,N)
         REAL*8, dimension(20) :: RPAR
         INTEGER, dimension(20) :: IPAR
         DO I=1,N
           AM(1,I)=1.0D0
         END DO
C        The first is also algebraic as is the integral form of the 
C        fractional differential equation         
         AM(1,1) = 0.0D0
         AM(1,2) = 0.0D0
         RETURN
         END
C        
        SUBROUTINE HMNEXP (T,ALPHA,EPSIL,DELTA,HEXP,MEXP,NEXP)
        IMPLICIT REAL*8 (A-H,O-Z)
C --- COMPUTES THE VALUES OF HEXP,MEXP,NEXP FOR THE FRACTIONAL DERIVATIVE
        MALPH=INT(ALPHA)
        ALPHA0=ALPHA-MALPH
        DELTA=(GAMMA(1+ALPHA)*EPSIL)**(1.D0/ALPHA)
        PI=4.D0*ATAN(1.D0)
        A=(PI/2.D0)*(1.D0+(1-ALPHA0)/((2-ALPHA0)*LOG(EPSIL)))
        HEXP=2.D0*PI*A/LOG(1.D0+(2.D0/EPSIL)*COS(A)**(ALPHA0-1))
        XL= (GAMMA(2-ALPHA0)*EPSIL)**(1.D0/(1-ALPHA0))
        XS= -LOG(GAMMA(1-ALPHA0)*EPSIL)
        NEXP = INT(LOG(XS/DELTA)/HEXP)+1
        MEXP = INT(LOG(XL/T)/HEXP)-1
        RETURN
        END
