C * * * * * * * * * * * * * * * * * * * * * * * * *
C --- DRIVER FOR RADAU5: FRACTIONAL LAPLACIAN  0 < alpha < 1
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 :: NDIM=61000,NDM=300
        PARAMETER (LWORK=NDIM*(11*NDM+41)+20,LIWORK=3*NDIM+20)               
        REAL*8, dimension(NDM+20) :: RPAR
        INTEGER, dimension(4) :: IPAR
        DIMENSION Y(NDIM),WORK(LWORK),IWORK(LIWORK),YEX(NDIM)
        DIMENSION ATOL(NDIM),RTOL(NDIM)
        EXTERNAL FGAMI,JGAMI,JGAMIEX,MAS,SOLOUT
        REAL*4 start,finish
C --- PARAMETERS OF THE PROBLEM, ALPHA BETWEEN 0 and 1
        ALPHA=1.D0/3.D0
        BETA =5.D0/3.D0
        RPAR(1)=ALPHA
        RPAR(2)=BETA
C --- DESIRED ACCURACY
        EPSIL=1.0D-6
        TOL=EPSIL
C --- END OF INTEGRATION
        T=1.D1
        ND=300
        ND=30
        CALL HMNEXP (T,ALPHA,EPSIL,DELTA,HEXP,MEXP,NEXP)
        write (6,*) '  '
        WRITE(*,*) 'DELTA,HEXP,MEXP,NEXP,N ', delta, HEXP, 
     1                       MEXP, NEXP, ND*(NEXP-MEXP+2)
        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
           IPAR(3)=ND
C          ND equations + ND integrals
           N = ND + ND + MI*ND*(NEXP-MEXP) 
           IF (ILIN.EQ.2) N = N + 2*ND
           IF (N.GT.NDIM.OR.ND.GT.NDM) THEN
              WRITE (6,*) ' NDIM or NDM too small',N,ND
              STOP
           END IF
C --- INITIAL VALUES
           X=0.0D0 
           XEND=T
           Y=0.0D0
           DX=1.D0/(ND+1)
C          INITIAL DATA           
           DO K=1,ND
               XK=K*DX
               Y(K)=0.5D0*XK*(1-XK)
               RPAR(K+10)=Y(K)
C              Integrals I_k
               Y(K+ND)=0.D0
           END DO
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=2*ND+1,N
               ATOL(I)=FACTOL
               RTOL(I)=FACTOL
             END DO
           END IF
           H=0.01D0
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
             MLJAC=2*ND
             MUJAC=0
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
C --- exact solution
           DO I=1,ND
               XI=I*DX
               YEX(I)=0.5D0*XI*(1-XI)*(T**BETA+1.D0)
           END DO
           ERR=0.0D0
           DO I=1,ND
             ERR = ERR + ((Y(I)-YEX(I))/YEX(I))**2
           END DO
           ERR=SQRT(ERR)
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(1000) :: RPAR
        INTEGER, dimension(20) :: IPAR
        DIMENSION Y(N),CONT(LRC)
        COMMON /INTERN/XOUT
        ND=IPAR(3)
        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)
        INTEGER, PARAMETER :: DP=kind(1D0) 
        REAL*8, dimension(1000) :: RPAR
        INTEGER, dimension(20) :: IPAR
        REAL(kind=DP), dimension(:), allocatable :: G,GAM,CI
        DIMENSION Y(N),F(N)
        ALPHA = RPAR(1)
        BETA  = RPAR(2)
        ETA   = BETA-ALPHA
        C = GAMMA(BETA+1)/GAMMA(BETA+1-ALPHA)
        PI=4.0D0*ATAN(1.0D0)
        MEXP=IPAR(1)
        NEXP=IPAR(2)
        LEXP=NEXP-MEXP
        HEXP=RPAR(6) 
        ND=IPAR(3)
        ALLOCATE(G(ND),GAM(LEXP),CI(LEXP)) 
        DX=1.D0/(ND+1)
        DX2=DX*DX
C       THE G FUNCTIONS        
        X1=DX
        G(1)=(-2.D0*Y(1)+Y(2))/DX2 
     1      + 0.5D0*X1*(1-X1)*C*X**ETA + X**BETA + 1.D0
        DO K=2,ND-1
            XK=K*DX
            G(K)=(Y(K-1)-2.D0*Y(K)+Y(K+1))/DX2 
     1          + 0.5D0*XK*(1-XK)*C*X**ETA + X**BETA + 1.D0
        END DO
        XND=ND*DX
        G(ND)=(Y(ND-1)-2.D0*Y(ND))/DX2 
     1      + 0.5D0*XND*(1-XND)*C*X**ETA + X**BETA + 1.D0
C       COMPUTATION OF COEFFICIENTS AND EXPONENTS        
        DO I=1,LEXP
           IND=MEXP+(I-1)
           GAM(I) =  EXP(IND*HEXP) 
           CI(I)  =  EXP((1-ALPHA)*IND*HEXP)
        END DO
C        
        II=2*ND+1
        DO K=1,ND
         APPINT=0.D0   
         DO I=1,LEXP
           F(II)   = - GAM(I)*Y(II) + G(K)
           APPINT =  APPINT + CI(I)*Y(II)
           II=II+1
         END DO
         APPINT=APPINT*HEXP*SIN(PI*ALPHA)/(PI)
         F(K) = RPAR(K+10) + Y(K+ND) - Y(K)
         F(K+ND) = Y(K+ND)-APPINT
        END DO 
        DEALLOCATE(G,GAM,CI)
        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)
        INTEGER, PARAMETER :: DP=kind(1D0) 
        REAL*8, dimension(1000) :: RPAR
        INTEGER, dimension(20) :: IPAR
        DIMENSION Y(N),DFY(LDFY,N)
        REAL(kind=DP), dimension(:), allocatable :: GAM,CI
        ALPHA = RPAR(1)
        PI=4.0D0*ATAN(1.0D0)
        MEXP=IPAR(1)
        NEXP=IPAR(2)
        LEXP=NEXP-MEXP
        HEXP=RPAR(6) 
        ND=IPAR(3)
        ALLOCATE(GAM(LEXP),CI(LEXP)) 
        DX=1.D0/(ND+1)
        DX2=DX*DX
C       COMPUTATION OF COEFFICIENTS AND EXPONENTS        
        DO I=1,LEXP
           IND=MEXP+(I-1)
           GAM(I) =  EXP(IND*HEXP) 
           CI(I)     =  EXP((1-ALPHA)*IND*HEXP)
        END DO
C
        DFY = 0.0D0     
        DO I=1,ND
            DFY(I,I) = -1.D0
            DFY(I,I+ND) = 1.D0
            DFY(I+ND,I+ND) = 1.D0
        END DO
        FAC=HEXP*SIN(PI*ALPHA)/(PI)
C
        II=2*ND+1
        K=1
        DO I=1,LEXP
          DFY(ND+K,II) = -CI(I)*FAC
          DFY(II,II) = - GAM(I)
          DFY(II,K) = -2.D0/DX2
          DFY(II,K+1) = 1.D0/DX2
          II=II+1
        END DO
        DO K=2,ND-1
         DO I=1,LEXP
          DFY(ND+K,II) = -CI(I)*FAC
          DFY(II,II) = - GAM(I)
          DFY(II,K-1) = 1.D0/DX2
          DFY(II,K) = -2.D0/DX2
          DFY(II,K+1) = 1.D0/DX2
          II=II+1
         END DO
        END DO
        K=ND
        DO I=1,LEXP
          DFY(ND+K,II) = -CI(I)*FAC
          DFY(II,II) = - GAM(I)
          DFY(II,K-1) = 1.D0/DX2
          DFY(II,K) = -2.D0/DX2
          II=II+1
        END DO
        DEALLOCATE(GAM,CI)
        RETURN
        END
C ==================================================================    
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)
        INTEGER, PARAMETER :: DP=kind(1D0) 
        REAL*8, dimension(1000) :: RPAR
        INTEGER, dimension(20) :: IPAR
        DIMENSION Y(N),DFY(LDFY,N)
        REAL(kind=DP), dimension(:), allocatable :: GAM,CI
        ALPHA = RPAR(1)
        PI=4.0D0*ATAN(1.0D0)
        MEXP=IPAR(1)
        NEXP=IPAR(2)
        LEXP=NEXP-MEXP
        HEXP=RPAR(6) 
        ND=IPAR(3)
        ALLOCATE(GAM(LEXP),CI(LEXP)) 
        DX=1.D0/(ND+1)
        DX2=DX*DX
C       COMPUTATION OF COEFFICIENTS AND EXPONENTS        
        DO I=1,LEXP
           IND=MEXP+(I-1)
           GAM(I) =  EXP(IND*HEXP) 
           CI(I)     =  EXP((1-ALPHA)*IND*HEXP)
        END DO
C
        DFY = 0.0D0     
        DO I=1,ND
            DFY(I,I) = -1.D0
            DFY(I,I+ND) = 1.D0
            DFY(I+ND,I+ND) = 1.D0
        END DO
        FAC=HEXP*SIN(PI*ALPHA)/(PI)
C
        II=2*ND+1
        K=1
C       First integral DF1/Dsum1        
        DFY(ND+K,II) =  -1.0D0
C       Dg1/Dy 
        II=II+1
        DFY(K,II) = -2.D0/DX2
        DFY(K+1,II) = 1.D0/DX2
        II=II+1
        DO I=1,LEXP
          DFY(1,II) = CI(I)*FAC
          DFY(2,II) = - GAM(I)
          II=II+1
        END DO
C       Separator        
        DFY(3,II-1)=-1.0D0
        DO K=2,ND-1
C        k-th integral DFk/Dsumk        
         DFY(ND+K,II) =  -1.0D0
C        Dgk/Dy 
         II=II+1
         DFY(K-1,II) =  1.D0/DX2
         DFY(K,  II) = -2.D0/DX2
         DFY(K+1,II) =  1.D0/DX2
         II=II+1
         DO I=1,LEXP
          DFY(1,II) = CI(I)*FAC
          DFY(2,II) = - GAM(I)
          II=II+1
         END DO
C       Separator        
         DFY(3,II-1)=-1.0D0
        END DO
        K=ND
C       Last integral DFNd/DsumNd        
        DFY(ND+K,II) =  -1.0D0
C       DgNd/Dy 
        II=II+1
        DFY(K-1,II) =  1.D0/DX2
        DFY(K  ,II) = -2.D0/DX2
        II=II+1
        DO I=1,LEXP
          DFY(1,II) = CI(I)*FAC
          DFY(2,II) = - GAM(I)
          II=II+1
        END DO
C       Separator        
        DFY(3,II-1)=-ND
        DEALLOCATE(GAM,CI)
        RETURN
        END
C
         SUBROUTINE MAS(N,AM,LMAS,RPAR,IPAR)
         IMPLICIT REAL*8 (A-H,O-Z)
         DIMENSION AM(LMAS,N)
         REAL*8, dimension(1000) :: RPAR
         INTEGER, dimension(20) :: IPAR
         ND=IPAR(3)
         DO I=1,N
           AM(1,I)=1.0D0
         END DO
         DO I=1,2*ND
           AM(1,I)=0.0D0
         END DO
         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
