C * * * * * * * * * * * * * * * * * * * * * * * * *
C --- DRIVER FOR RADAU5: FRACTIONAL BRUSSELATOR PROBLEM
C     version of 27 May 2025
C * * * * * * * * * * * * * * * * * * * * * * * * *
        include '../../radau5.f'
        include '../../decsol.f'
        include 'dc_sumexp.f'
        IMPLICIT REAL*8 (A-H,O-Z)
        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
        T=220.0D0
        alpha1=1.3D0
        alpha2=0.8D0
        RPAR(1)=alpha1
        RPAR(2)=alpha2
        A=1.0D0
        B=3.0D0
        RPAR(11)=A
        RPAR(12)=B
C --- DESIRED ACCURACY
        EPSIL=1.0D-8
        TOL=EPSIL
        ALPH01=ALPHA1-1.0D0
        CALL HMNEXP (T,ALPH01,EPSIL,DELTA1,HEXP1,MEXP1,NEXP1)
        write (6,*) '  '
        WRITE(*,*) 'DELTA1,HEXP1,MEXP1,NEXP1 ', delta1, HEXP1, 
     1                                          MEXP1, NEXP1
        RPAR(6)=HEXP1
        IPAR(1)=MEXP1
        IPAR(2)=NEXP1
        CALL HMNEXP (T,ALPHA2,EPSIL,DELTA2,HEXP2,MEXP2,NEXP2)
        write (6,*) '  '
        WRITE(*,*) 'DELTA2,HEXP2,MEXP2,NEXP2 ', delta2, HEXP2, 
     1                                          MEXP2, NEXP2
        RPAR(7)=HEXP2
        IPAR(3)=MEXP2
        IPAR(4)=NEXP2
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 = 2,2
           write (6,*) '  '
           MI=1
           ND=2
C          ND = 2 equations
           N = ND + MI*(NEXP1-MEXP1) + MI*(NEXP2-MEXP2) 
           IF (ILIN.EQ.2) N = N + 2 + 2
C --- INITIAL VALUES
           X=0.0D0 
           XEND=T
           Y=0.0D0
           Y(1)=1.2D0
           Y(2)=2.8D0
           DY01=1.0D0
           RPAR(9)=DY01
           RPAR(10)=Y(2)
C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM? NO
           IMAS=1
           MLMAS=0
           MUMAS=0
C --- OUTPUT ROUTINE IS USED DURING INTEGRATION?    
           IOUT=1
           IOUT=0
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
             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) = 1.0097684171D+00
           YEX(2) = 2.1581264031D+00
           ERR=0.0D0
           DO I=1,2
             ERR = ERR + ((Y(I)-YEX(I))/YEX(I))**2
           END DO
           ERR=SQRT(ERR/2)
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:2),ERR
 99        FORMAT(1X,'X =',F10.3,'    Y =',2E20.12,'    ERR = ',E9.3)
           WRITE (6,91) (IWORK(J),J=14,20)
 91        FORMAT(' fcn=',I6,' jac=',I6,' step=',I6,' accpt=',I6,
     &        ' rejct=',I4,' dec=',I6,' sol=',I6)
        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 (9,999) XOUT,CONTR5(1,XOUT,CONT,LRC),
     1         CONTR5(2,XOUT,CONT,LRC),CONTR5(3,XOUT,CONT,LRC),
     2         CONTR5(4,XOUT,CONT,LRC), NR-1        
              XOUT=XOUT+1.0D-2
              GOTO 10
           END IF
        END IF
 99     FORMAT(1X,'X =',F5.1,'    Y =',E18.10,'    NSTEP = ',I4)
 999    FORMAT(1X,'    ',F7.3,'     ',4E18.10,'     ',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        
        DY01=RPAR(9)
        Y02=RPAR(10)
        ALPH1 = RPAR(1)
        ALPH2 = RPAR(2)
        K1 = INT(ALPH1)
        K2 = INT(ALPH2)
        PI=4.0D0*ATAN(1.0D0)
        MEXP1=IPAR(1)
        NEXP1=IPAR(2)
        LEXP1=NEXP1-MEXP1
        HEXP1=RPAR(6) 
        MEXP2=IPAR(3)
        NEXP2=IPAR(4)
        LEXP2=NEXP2-MEXP2
        HEXP2=RPAR(7)
        A=RPAR(11)
        B=RPAR(12)
C        
        F1 = A - (B+1)*Y(1) + Y(1)*Y(1)*Y(2)
        F2 = B*Y(1) - Y(1)*Y(1)*Y(2)
C        
        ALPH0=ALPH1-K1
        APPINT1=0.0D0
        II=3
        DO I=1,LEXP1
           IND=MEXP1+(I-1)
           GAMMAI =  EXP(IND*HEXP1) 
           CI     =  EXP((1-ALPH0)*IND*HEXP1)
           F(II)   = - GAMMAI*Y(II) + F1
           APPINT1 =  APPINT1 + CI*Y(II)
           II=II+1
        END DO
        APPINT1=APPINT1*HEXP1*SIN(PI*ALPH0)/PI
C        
        ALPH0=ALPH2-K2
        APPINT2=0.0D0
        DO I=1,LEXP2    
           IND=MEXP2+(I-1)
           GAMMAI =  EXP(IND*HEXP2) 
           CI     =  EXP((1-ALPH0)*IND*HEXP2)
           F(II)   = - GAMMAI*Y(II) + F2
           APPINT2 =  APPINT2 + CI*Y(II)
           II=II+1
        END DO
        APPINT2=APPINT2*HEXP2*SIN(PI*ALPH0)/PI           
C
        F(1) = DY01 + APPINT1 
        F(2) = Y02 + APPINT2 - 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)
        PI=4.0D0*ATAN(1.0D0)
        alph1 = RPAR(1)
        alph2 = RPAR(2)
        k1 = INT(alph1)
        k2 = INT(alph2)
        MEXP1=IPAR(1)
        NEXP1=IPAR(2)
        LEXP1=NEXP1-MEXP1
        HEXP1=RPAR(6) 
        MEXP2=IPAR(3)
        NEXP2=IPAR(4)
        LEXP2=NEXP2-MEXP2
        HEXP2=RPAR(7)   
        A=RPAR(11)
        B=RPAR(12)
C        
        DFY = 0.0D0    
        DFY(2,2) = -1.D0        
        ALPH0=ALPH1-K1
        FAC1=HEXP1*SIN(PI*ALPH0)/PI
        II=3
        DO I=1,LEXP1
          IND=MEXP1+(I-1)
          GAMMAI =  EXP(IND*HEXP1) 
          CI =      EXP((1-alph0)*IND*HEXP1)
          DFY(1,II) = CI*FAC1
          DFY(II,II) = - GAMMAI
          DFY(II,1) = -(B+1) + 2.D0*Y(1)*Y(2)
          DFY(II,2) =  Y(1)*Y(1)
          II=II+1
        END DO
        ALPH0=ALPH2-K2
        FAC2=HEXP2*SIN(PI*ALPH0)/PI
        DO I=1,LEXP2    
           IND=MEXP2+(I-1)
           GAMMAI =  EXP(IND*HEXP2) 
           CI     =  EXP((1-alph0)*IND*HEXP2)
           DFY(2,II) = CI*FAC2
           DFY(II,II) = - GAMMAI
           DFY(II,1) = B - 2.D0*Y(1)*Y(2)
           DFY(II,2) = -Y(1)*Y(1)
           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)
        PI=4.0D0*ATAN(1.0D0)
        alph1 = RPAR(1)
        alph2 = RPAR(2)
        k1 = INT(alph1)
        k2 = INT(alph2)
        MEXP1=IPAR(1)
        NEXP1=IPAR(2)
        LEXP1=NEXP1-MEXP1
        HEXP1=RPAR(6) 
        MEXP2=IPAR(3)
        NEXP2=IPAR(4)
        LEXP2=NEXP2-MEXP2
        HEXP2=RPAR(7)   
        A=RPAR(11)
        B=RPAR(12)
C        
        DFY = 0.0D0    
        DFY(2,2) = -1.D0  
C       First integral DF1/Dsum1  
        II=3      
        DFY(1,II) = 1.0D0
        II=II+1
C       Dg1/Dy 
        DFY(1,II) = -(B+1) + 2.D0*Y(1)*Y(2)
        DFY(2,II) = Y(1)*Y(1)
        II=II+1
        ALPH0=ALPH1-K1
        FAC1=HEXP1*SIN(PI*ALPH0)/PI
        DO I=1,LEXP1
          IND=MEXP1+(I-1)
          GAMMAI =  EXP(IND*HEXP1) 
          CI =      EXP((1-alph0)*IND*HEXP1)
          DFY(1,II) = CI*FAC1
          DFY(2,II) = - GAMMAI
          II=II+1
        END DO
C       Separator        
        DFY(3,II-1)=-1.0D0
C       Second integral DF/Dsum2
        DFY(2,II) = 1.D0
        II = II+1
C       Dg2/Dy        
        DFY(1,II) = B - 2.D0*Y(1)*Y(2)
        DFY(2,II) = -Y(1)*Y(1)        
        II=II+1
        ALPH0=ALPH2-K2
        FAC2=HEXP2*SIN(PI*ALPH0)/PI
        DO I=1,LEXP2
          IND=MEXP2+(I-1)
          GAMMAI =  EXP(IND*HEXP2) 
          CI =      EXP((1-alph0)*IND*HEXP2)
          DFY(1,II) = CI*FAC2
          DFY(2,II) = - GAMMAI
          II=II+1
        END DO
C       Separator        
        DFY(3,II-1)=-2.0D0
        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
         AM(1,1) = 1.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
