C * * * * * * * * * * * * * * * * * * * * * * * * *
C --- DRIVER FOR RADAU5: FRACTIONAL SYSTEM OF LAPLACIANS  7-diagonal ordering
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=5000000,NDM=20000
        PARAMETER (LWORK=NDIM*(28)+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
        open(unit=10, file="exact.txt")
C --- PARAMETERS OF THE PROBLEM
        ALPHA=1.D0/2.D0
        AK1=1.0D0
        AK2=2.0D0
        AK3=3.0D0
        AK=0.5D0
        RPAR(1)=ALPHA
        RPAR(2)=AK1
        RPAR(3)=AK2
        RPAR(4)=AK3
        RPAR(5)=AK
C --- DESIRED ACCURACY
        EPSIL=1.0D-6
        TOL=EPSIL
C --- END OF INTEGRATION
        T=30.0D0
        ND=1000
c        ND=7
        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 --- FULL LINEAR ALGEBRA (ILIN = 1) AND ADAPTED TO SUMEXP (ILIN = 2)
        DO ILIN = 2,2
           write (6,*) '  '
           MI=1
           IPAR(3)=ND
           ND3=3*ND
C          ND equations 
           N = ND3 + MI*ND3*(NEXP-MEXP) 
           IF (ILIN.GE.2) N = N + 2*ND3
           IF (N.GT.NDIM.OR.ND.GT.NDM) THEN
              WRITE (6,*) ' NDIM or NDM too small',N,ND
c              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
            K3=3*K
               XK=K*DX
               Y(K3-2)=0.5D0*XK*(1-XK)
               Y(K3-1)=1.0D0*XK**2*(1-XK)
               Y(K3)=1.5D0*XK*(1-XK)**2
           END DO
           DO K=1,ND3
               RPAR(K+10)=Y(K)
           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
           H=0.001D0
C --- SET DEFAULT VALUES 
           DO I=1,20
             IWORK(I)=0
             WORK(I)=0.D0
           END DO
c           IWORK(2)=3  ! maximal number of steps
C ==================================================================
C          FULL LIN. ALGEBRA OPTION           
C ==================================================================           
           IF (ILIN.EQ.1) THEN           
             IJAC=1
c             IJAC=0
             MUJAC=0
             MLJAC=N
             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 (BANDED IN DIMENSION ND))
C ==================================================================           
           IF (ILIN.EQ.2) THEN
             IJAC=1
             MUJAC=3
             MLJAC=3
             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 ==================================================================           
c
           do i=1,nd3
             read (10,*) ii,yex(i)
           end do
           ERR=0.0D0
           DO I=1,ND3
             ERR = ERR + ((Y(I)-YEX(I))/YEX(I))**2
           END DO
           ERR=SQRT(ERR/ND3)
C --- PRINT FINAL SOLUTION
           write (6,77) ILIN,TOL,finish-start
 77        FORMAT(1X,'radau5, lin-alg = ',I1,'   TOL = ',E9.3,
     &           '    time = ',E9.3)
           WRITE (6,99) X,Y(7),Y(8),Y(9),ERR
 99        FORMAT(1X,'X =',F10.3,'    Y =',3E20.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
        close (10)
        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)
        AK1 = RPAR(2)
        AK2 = RPAR(3)
        AK3 = RPAR(4)
        AK = RPAR(5)
        PI=4.0D0*ATAN(1.0D0)
        MEXP=IPAR(1)
        NEXP=IPAR(2)
        LEXP=NEXP-MEXP
        HEXP=RPAR(6) 
        ND=IPAR(3)
        ND2=2*ND
        ND3=3*ND
        ALLOCATE(G(ND3),GAM(LEXP),CI(LEXP)) 
        DX=1.D0/(ND+1)
        DX2=DX*DX/AK
C       THE G FUNCTIONS        
        G(1)=(-2.D0*Y(1)+Y(4))/DX2 
     1      - AK1*Y(1)*Y(2) + (AK2+AK3)*Y(3)
        G(2)=(-2.D0*Y(2)+Y(5))/DX2 
     1      - AK1*Y(1)*Y(2) + AK2*Y(3)
        G(3)=(-2.D0*Y(3)+Y(6))/DX2 
     1      + AK1*Y(1)*Y(2) - (AK2+AK3)*Y(3)
        DO K=2,ND-1
            K3=3*K
            G(K3-2)=(Y(K3-5)-2.D0*Y(K3-2)+Y(K3+1))/DX2 
     1          - AK1*Y(K3-2)*Y(K3-1) + (AK2+AK3)*Y(K3)
            G(K3-1)=(Y(K3-4)-2.D0*Y(K3-1)+Y(K3+2))/DX2 
     1          - AK1*Y(K3-2)*Y(K3-1) + AK2*Y(K3)
            G(K3)=(Y(K3-3)-2.D0*Y(K3)+Y(K3+3))/DX2 
     1          + AK1*Y(K3-2)*Y(K3-1) - (AK2+AK3)*Y(K3)
        END DO
        G(ND3-2)=(Y(ND3-5)-2.D0*Y(ND3-2))/DX2 
     1      - AK1*Y(ND3-2)*Y(ND3-1) + (AK2+AK3)*Y(ND3)
        G(ND3-1)=(Y(ND3-4)-2.D0*Y(ND3-1))/DX2 
     1      - AK1*Y(ND3-2)*Y(ND3-1) + AK2*Y(ND3)
        G(ND3)=(Y(ND3-3)-2.D0*Y(ND3))/DX2 
     1      + AK1*Y(ND3-2)*Y(ND3-1) - (AK2+AK3)*Y(ND3)
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
        FAC = HEXP*SIN(PI*ALPHA)/(PI)
C        
        II=ND3+1
        DO KK=1,ND
         DO ILAP=1,3
          K=3*(KK-1)+ILAP
          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*FAC
          F(K) = RPAR(K+10) + APPINT - Y(K)
         END DO 
        END DO
        DO I=II,N
          F(I)=0.0D0
        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)
        AK1 = RPAR(2)
        AK2 = RPAR(3)
        AK3 = RPAR(4)
        AK = RPAR(5)
        PI=4.0D0*ATAN(1.0D0)
        MEXP=IPAR(1)
        NEXP=IPAR(2)
        LEXP=NEXP-MEXP
        HEXP=RPAR(6) 
        ND=IPAR(3)
        ND2=2*ND
        ND3=3*ND
        ALLOCATE(GAM(LEXP),CI(LEXP)) 
        DX=1.D0/(ND+1)
        DX2=DX*DX/AK
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,ND3
            DFY(I,I) = -1.D0
        END DO
C
        FAC=HEXP*SIN(PI*ALPHA)/(PI)
        II=ND3+1
        DO KK=1,ND
         K=3*(KK-1)+1
         DO I=1,LEXP
          DFY(K,II) = CI(I)*FAC
          DFY(II,II) = - GAM(I)
            IF (KK.GT.1) DFY(II,K-3) = 1.D0/DX2
            DFY(II,K) = -2.D0/DX2 - AK1*Y(K+1)
            DFY(II,K+1) = - AK1*Y(K)
            DFY(II,K+2) = (AK2+AK3)
            IF (KK.LT.ND) DFY(II,K+3) = 1.D0/DX2
            II=II+1
         END DO
         DO I=1,LEXP
          DFY(K+1,II) = CI(I)*FAC
          DFY(II,II) = - GAM(I)
            IF (KK.GT.1) DFY(II,K-2) = 1.D0/DX2
            DFY(II,K) = - AK1*Y(K+1)
            DFY(II,K+1) = -2.D0/DX2 - AK1*Y(K)
            DFY(II,K+2) = AK2
            IF (KK.LT.ND) DFY(II,K+4) = 1.D0/DX2
            II=II+1
         END DO
         DO I=1,LEXP
          DFY(K+2,II) = CI(I)*FAC
          DFY(II,II) = - GAM(I)
            IF (KK.GT.1) DFY(II,K-1) = 1.D0/DX2
            DFY(II,K) = AK1*Y(K+1)
            DFY(II,K+1) = AK1*Y(K)
            DFY(II,K+2) = -2.D0/DX2 - (AK2+AK3)
            IF (KK.LT.ND) DFY(II,K+5) = 1.D0/DX2
          II=II+1
         END DO
        END DO
        DEALLOCATE(GAM,CI)
        RETURN
        END
C
C ==================================================================
C          SUMEXP LIN. ALGEBRA OPTION  (BANDED IN DIMENSION ND)
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)
        AK1 = RPAR(2)
        AK2 = RPAR(3)
        AK3 = RPAR(4)
        AK = RPAR(5)
        PI=4.0D0*ATAN(1.0D0)
        MEXP=IPAR(1)
        NEXP=IPAR(2)
        LEXP=NEXP-MEXP
        HEXP=RPAR(6) 
        ND=IPAR(3)
        ND2=2*ND
        ND3=3*ND
        ALLOCATE(GAM(LEXP),CI(LEXP)) 
        DX=1.D0/(ND+1)
        DX2=DX*DX/AK
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
        MLJAC = 3
        MUJAC = 3
        MDIAG = MUJAC + 1
        DFY = 0.0D0
C      number of integrals
        DFY(1,1) = ND3
        DO I=1,ND3
            DFY(MDIAG,I) = -1.D0
        END DO
        FAC=HEXP*SIN(PI*ALPHA)/(PI)
        II=ND3+1
C
        DO KK=1,ND
         K=3*(KK-1)+1
         DFY(MDIAG,II) =  1.0D0
            II=II+1
            IF (KK.GT.1) DFY(MDIAG-3,II) = 1.D0/DX2
            DFY(MDIAG,II) = -2.D0/DX2 - AK1*Y(K+1)
            DFY(MDIAG+1,II) = - AK1*Y(K)
            DFY(MDIAG+2,II) = (AK2+AK3)
            IF (KK.LT.ND) DFY(MDIAG+3,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)=-ND3
c
         DFY(MDIAG,II) =  1.0D0
            II=II+1
            IF (KK.GT.1) DFY(MDIAG-3,II) = 1.D0/DX2
            DFY(MDIAG-1,II) = - AK1*Y(K+1)
            DFY(MDIAG,II) = -2.D0/DX2 - AK1*Y(K)
            DFY(MDIAG+1,II) = AK2
            IF (KK.LT.ND) DFY(MDIAG+3,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)=-ND3
c
         DFY(MDIAG,II) =  1.0D0
            II=II+1
            IF (KK.GT.1) DFY(MDIAG-3,II) = 1.D0/DX2
            DFY(MDIAG-2,II) = AK1*Y(K+1)
            DFY(MDIAG-1,II) = AK1*Y(K)
            DFY(MDIAG,II) = -2.D0/DX2 - (AK2+AK3)
            IF (KK.LT.ND) DFY(MDIAG+3,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)=-ND3
        END DO
        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,3*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
