PROGRAM TEST
    USE M_PENTA
    CALL EX05
CONTAINS
!
    SUBROUTINE EX01    
    INTEGER, PARAMETER :: NN = 10
    INTEGER :: N, IER
    REAL(8) :: A(NN), B(NN), C(NN), D(NN), E(NN), F(NN), R(NN)
    REAL(8) :: A0(NN), B0(NN), C0(NN), D0(NN), E0(NN), F0(NN)
    N = NN
    A = DBLE((/0,0,-2,-5,1,2,3,1,3,4/))  
    B = DBLE((/4,2,2,1,7,-3,0,2,2,-3/))
    C = DBLE((/1,3,1,2,-5,4,7,3,-5,2/))
    D = DBLE((/2,4,3,-1,2,-1,-8,-4,1,-2/))
    E = DBLE((/3,5,2,3,-2,2,1,2,0,0/))
    F = DBLE((/10,14,6,0,3,4,3,4,1,1/)) 
    A0 = A
    B0 = B
    C0 = C
    D0 = D
    E0 = E
    F0 = F 
    CALL PENTA( N,A,B,C,D,E,F,IER)
    PRINT *,'IER=',IER
    PRINT *,F
    CALL PENRES1( N, A0, B0, C0, D0, E0, F0, F, R)
    PRINT *,R
    END SUBROUTINE EX01
!
    SUBROUTINE EX02a   
        INTEGER, PARAMETER :: MM = 8
        INTEGER, PARAMETER :: NN = 100000
        INTEGER :: M, N, IER, K
        REAL(8) :: A(MM,MM,NN), B(MM,MM,NN), C(MM,MM,NN), D(MM,MM,NN), &
        &   E(MM,MM,NN), F(MM,NN), R(MM,NN)
        REAL(8) :: A0(MM,MM,NN), B0(MM,MM,NN), C0(MM,MM,NN), D0(MM,MM,NN), &
        &   E0(MM,MM,NN), F0(MM,NN), RMAX,RAVE,BETA
        real t1, t2
        BETA=4.0d0
        M = MM
        N = NN
        A = 0.0D0
        DO K = 1, M
           a(K,K,1:NN) = 1.0d0 !DBLE(2*MM) + C(K,K,1:NN)
        ENDDO        
        B=2*A
        C=A
        D=-A
        E=-2*A   
        DO K = 1, M
            C(K,K,1:NN) = DBLE(BETA) + C(K,K,1:NN)
        ENDDO    
        A0 = A
        B0 = B
        C0 = C
        D0 = D
        E0 = E
        F0 = 1.0D0
        R  = 0.0D0 
        
        CALL PENRES2(M, N, A0, B0, C0, D0, E0, R, F0, F)
        F0 = F
      !  PRINT *,F
        CALL CPU_TIME(T1)
        CALL PENTA( M,N,A,B,C,D,E,F,IER) !,(/1.0D0,1.0D0,1.0D0,1.0D0/))
        CALL CPU_TIME(T2)
        PRINT *,'IER=',IER,' SEC:',T2-T1
       ! PRINT *,'F=',F
        CALL PENRES2(M, N, A0, B0, C0, D0, E0, F0, F, R)
      !  PRINT *,'R=',R
         RMAX = MAXVAL(ABS(R))
         RAVE = SUM(ABS(R))/(N*DBLE(M))
         PRINT *,'RMAX=',RMAX,' RAVE=',RAVE
         F0=1.0d0
         PRINT *,'ERR=',MAXVAL(ABS(F-F0)) 
    END SUBROUTINE EX02a
    SUBROUTINE EX02   
        INTEGER, PARAMETER :: MM = 8
        INTEGER, PARAMETER :: NN = 100000
        INTEGER :: M, N, IER, K
        REAL(8) :: A(MM,MM,NN), B(MM,MM,NN), C(MM,MM,NN), D(MM,MM,NN), &
        &   E(MM,MM,NN), F(MM,NN), R(MM,NN)
        REAL(8) :: A0(MM,MM,NN), B0(MM,MM,NN), C0(MM,MM,NN), D0(MM,MM,NN), &
        &   E0(MM,MM,NN), F0(MM,NN), RMAX,RAVE
        real t1, t2
        M = MM
        N = NN
        CALL RANDOM_NUMBER(A)  
        CALL RANDOM_NUMBER(B)  
        CALL RANDOM_NUMBER(C)  
        CALL RANDOM_NUMBER(D)  
        CALL RANDOM_NUMBER(E) 
      !  A = 0.0D0
      
        DO K = 1, M
        !    a(K,K,1:NN) = 1.0d0 !DBLE(2*MM) + C(K,K,1:NN)
        ENDDO           
        DO K = 1, M
           ! a(K,K,1:NN) = DBLE(1) + a(K,K,1:NN)
           ! b(K,K,1:NN) = DBLE(1) + b(K,K,1:NN)
           ! d(K,K,1:NN) = DBLE(1) + d(K,K,1:NN)
           ! e(K,K,1:NN) = DBLE(1) + e(K,K,1:NN)

            C(K,K,1:NN) = DBLE(4*M) + C(K,K,1:NN)
        ENDDO    
        A0 = A
        B0 = B
        C0 = C
        D0 = D
        E0 = E
        F0 = 1.0D0
        R  = 0.0D0 
        
        CALL PENRES2(M, N, A0, B0, C0, D0, E0, R, F0, F)
        F0 = F
      !  PRINT *,F
        CALL CPU_TIME(T1)
        CALL PENTA( M,N,A,B,C,D,E,F,IER) !,(/1.0D0,1.0D0,1.0D0,1.0D0/))
        CALL CPU_TIME(T2)
        PRINT *,'IER=',IER,' SEC:',T2-T1
       ! PRINT *,'F=',F
        CALL PENRES2(M, N, A0, B0, C0, D0, E0, F0, F, R)
      !  PRINT *,'R=',R
         RMAX = MAXVAL(ABS(R))
         RAVE = SUM(ABS(R))/(N*DBLE(M))
         PRINT *,'RMAX=',RMAX,' RAVE=',RAVE
         F0=1.0d0
         PRINT *,'ERR=',MAXVAL(ABS(F-F0)) 
    END SUBROUTINE EX02 
!
!
    SUBROUTINE EX03   
        INTEGER, PARAMETER :: MM = 2
        INTEGER, PARAMETER :: NN = 2560
        INTEGER :: M, N, IER, K
        REAL(8) :: A(MM,MM,NN), B(MM,MM,NN), C(MM,MM,NN), D(MM,MM,NN), &
        &   E(MM,MM,NN), F(MM,NN), R(MM,NN)
        REAL(8) :: A0(MM,MM,NN), B0(MM,MM,NN), C0(MM,MM,NN), D0(MM,MM,NN), &
        &   E0(MM,MM,NN), F0(MM,NN), RMAX,RAVE
        REAL(8) :: H, PI, X
        PI = 3.1415926535897932384626433832795_8
        M = MM
        N = NN
        H = 1.0D0/DBLE(NN)
        A = 0.0D0
        B = A
        C = A
        D = A
        E = A        
        DO K = 1, M
            A(K,K,1:NN) = -1.0d0 
            B(K,K,1:NN) = 16.0d0 
            C(K,K,1:NN) = -30.0d0 
            D(K,K,1:NN) = 16.0d0 
            E(K,K,1:NN) = -1.0d0             
        ENDDO   
        C(1,2,1:NN) = 12.0D0*H**2  
        C(2,1,1:NN) = 12.0D0*H**2     
        DO K = 1, NN
            X = 2*PI*H*(K-1)
            F(1,K) = 12.0D0*H**2*(COS(X) - 4*PI**2*SIN(X)) 
            F(2,K) = 12.0D0*H**2*(SIN(X) - 4*PI**2*COS(X))
        ENDDO                      
        A0 = A
        B0 = B
        C0 = C
        D0 = D
        E0 = E
        F0 = F
        CALL PENTA( M,N,A,B,C,D,E,F,IER) !,(/1.0D0,1.0D0,1.0D0,1.D0/))
        PRINT *,'IER=',IER
        PRINT *,'F=',F
        OPEN(10,FILE='TEST.TXT',STATUS='UNKNOWN')
        DO K = 1, NN
            X = H*(K-1)        
            WRITE(10,'(13F16.8)') X,F(1,K),F(2,K),F(1,K)-SIN(2*PI*X),F(2,K)-COS(2*PI*X)
        ENDDO
        CLOSE(10)    
        CALL PENRES2(M, N, A0, B0, C0, D0, E0, F0, F, R)
        PRINT *,'R=',R
        RMAX = MAXVAL(ABS(R))
        RAVE = SUM(ABS(R))/(N*DBLE(M))
        PRINT *,'RMAX=',RMAX,' RAVE=',RAVE
        PAUSE
        DO K = 1, NN
            X = 2*PI*H*(K-1)
            R(1,K) = F(1,K) - SIN(X) 
            R(2,K) = F(2,K) - COS(X)
        ENDDO    
        PRINT *,'R=',R
        RMAX = MAXVAL(ABS(R))
        RAVE = SUM(ABS(R))/(N*DBLE(M))
        PRINT *,'RMAX=',RMAX,' RAVE=',RAVE
    END SUBROUTINE EX03 
!
!
    SUBROUTINE EX04
        INTEGER, PARAMETER :: MM = 2
        INTEGER, PARAMETER :: NN = 5
        INTEGER :: M, N, IER, K
        REAL(8) :: A(MM,MM,NN), B(MM,MM,NN), C(MM,MM,NN), D(MM,MM,NN), &
        &   E(MM,MM,NN), F(MM,NN), R(MM,NN)
        REAL(8) :: A0(MM,MM,NN), B0(MM,MM,NN), C0(MM,MM,NN), D0(MM,MM,NN), &
        &   E0(MM,MM,NN), F0(MM,NN), RMAX,RAVE
        M = MM
        N = NN

        A = 0.0D0
        C = 0.0D0
        DO K = 1, M
            A(K,K,1:NN) = 1.0d0 !DBLE(2*MM) + C(K,K,1:NN)
            C(K,K,1:NN) = 5.0d0 !DBLE(2*MM) + C(K,K,1:NN)            
        ENDDO           
        B = A
        D = A
        E = A
        A(:,:,1) = 0.0d0
        A(:,:,2) = 0.0d0
        E(:,:,N-1) = 0.0d0
        E(:,:,N) = 0.0d0                  
        A0 = A
        B0 = B
        C0 = C
        D0 = D
        E0 = E
        F0 = 1.0D0
        R  = 0.0D0 
        CALL PENRES2(M, N, A0, B0, C0, D0, E0, R, F0, F)
        F0 = F
        PRINT *,F
        CALL PENTA( M,N,A,B,C,D,E,F,IER) !,(/1.0D0,1.0D0,1.0D0,1.D0/))
        PRINT *,'IER=',IER
        PRINT *,'F=',F
        CALL PENRES2(M, N, A0, B0, C0, D0, E0, F0, F, R)
        PRINT *,'R=',R
         RMAX = MAXVAL(ABS(R))
         RAVE = SUM(ABS(R))/(N*DBLE(M))
         PRINT *,'RMAX=',RMAX,' RAVE=',RAVE
    END SUBROUTINE EX04    
!
!
    SUBROUTINE EX05
        INTEGER, PARAMETER :: MM = 2
        INTEGER, PARAMETER :: NN = 5
        INTEGER :: M, N, IER, K
        REAL(8) :: A(MM,MM,NN), B(MM,MM,NN), C(MM,MM,NN), D(MM,MM,NN), &
        &   E(MM,MM,NN), F(MM,NN), R(MM,NN)
        REAL(8) :: A0(MM,MM,NN), B0(MM,MM,NN), C0(MM,MM,NN), D0(MM,MM,NN), &
        &   E0(MM,MM,NN), F0(MM,NN), RMAX,RAVE
        M = MM
        N = NN

        A = 0.0D0
        C = 0.0D0
        A(:,:,1) = RESHAPE(DBLE((/1,1,1,-1/)),(/2,2/))
        B(:,:,1) = RESHAPE(DBLE((/-1,1,1,1/)),(/2,2/))
        C(:,:,1) = RESHAPE(DBLE((/5,1,1,5/)),(/2,2/))
        D(:,:,1) = RESHAPE(DBLE((/1,1,-1,1/)),(/2,2/))
        E(:,:,1) = RESHAPE(DBLE((/1,-1,1,1/)),(/2,2/))        
        DO K = 2, N
            A(:,:,K) = A(:,:,1)
            B(:,:,K) = B(:,:,1)
            C(:,:,K) = C(:,:,1)
            D(:,:,K) = D(:,:,1)
            E(:,:,K) = E(:,:,1)
        ENDDO    
        A(:,:,1) = 0.0d0
        A(:,:,2) = 0.0d0
        E(:,:,N-1) = 0.0d0
        E(:,:,N) = 0.0d0   
       ! PRINT *,D(1,1,3),D(1,2,3)
       ! PRINT *,D(2,1,3),D(2,2,3)
        
       ! PAUSE               
        A0 = A
        B0 = B
        C0 = C
        D0 = D
        E0 = E
        F0 = 1.0D0
        R  = 0.0D0 
        CALL PENRES2(M, N, A0, B0, C0, D0, E0, R, F0, F)
        F0 = F
        PRINT *,F
        CALL PENTA( M,N,A,B,C,D,E,F,IER,(/2.0D0,2.0D0,2.0D0,2.D0/))
        PRINT *,'IER=',IER
        PRINT *,'F=',F
        CALL PENRES2(M, N, A0, B0, C0, D0, E0, F0, F, R)
        PRINT *,'R=',R
         RMAX = MAXVAL(ABS(R))
         RAVE = SUM(ABS(R))/(N*DBLE(M))
         PRINT *,'RMAX=',RMAX,' RAVE=',RAVE
    END SUBROUTINE EX05
!    
!
    SUBROUTINE EX06
        INTEGER, PARAMETER :: MM = 7
        INTEGER, PARAMETER :: NN = 500
        INTEGER :: M, N, IER, K
        REAL(8) :: A(MM,MM,NN), B(MM,MM,NN), C(MM,MM,NN), D(MM,MM,NN), &
        &   E(MM,MM,NN), F(MM,NN), R(MM,NN)
        REAL(8) :: A0(MM,MM,NN), B0(MM,MM,NN), C0(MM,MM,NN), D0(MM,MM,NN), &
        &   E0(MM,MM,NN), F0(MM,NN), RMAX,RAVE
        REAL(4) :: T1, T2
        M = MM
        N = NN

        A = 0.0D0
        B = 1.8d0
        C = 1.0D0
        D = B
        E = A
       ! A(:,:,1) = RESHAPE(DBLE((/1,1,1,-1/)),(/2,2/))
       ! B(:,:,1) = RESHAPE(DBLE((/-1,1,1,1/)),(/2,2/))
       ! C(:,:,1) = RESHAPE(DBLE((/5,1,1,5/)),(/2,2/))
       ! D(:,:,1) = RESHAPE(DBLE((/1,1,-1,1/)),(/2,2/))
       ! E(:,:,1) = RESHAPE(DBLE((/1,-1,1,1/)),(/2,2/))        
        DO K = 1, M
            A(K,K,1:N) =  1.0D0
            B(K,K,1:N) = -7.2D0
        ENDDO    
        E = A
        D = B
        DO K = 1, M
        ENDDO    
        DO K = 1, M
            C(K,K,1:N) = 22.0D0
            IF (K < M ) THEN
                C(K,K+1,1:N)=-8.0D0
            ELSE
                C(K,1,1:N) = -8.0D0
            ENDIF
            IF (K == 1) THEN
                C(K,M,1:N) = -8.0D0
            ELSE
                C(K,K-1,1:N) = -8.0D0
            ENDIF                
        ENDDO    
        A0 = A
        B0 = B
        C0 = C
        D0 = D
        E0 = E
        F0 = 1.0D0
        R  = 0.0D0 
        
        CALL PENRES2(M, N, A0, B0, C0, D0, E0, R, F0, F)
        F0 = F
      !  PRINT *,F
        CALL CPU_TIME(T1)
        CALL PENTA( M,N,A,B,C,D,E,F,IER) !,(/1.0D0,1.0D0,1.0D0,1.0D0/))
        CALL CPU_TIME(T2)
        PRINT *,'IER=',IER,' SEC:',T2-T1
       ! PRINT *,'F=',F
        CALL PENRES2(M, N, A0, B0, C0, D0, E0, F0, F, R)
      !  PRINT *,'R=',R
         RMAX = MAXVAL(ABS(R))
         RAVE = SUM(ABS(R))/(N*DBLE(M))
         PRINT *,'RMAX=',RMAX,' RAVE=',RAVE
         F0=1.0d0
         PRINT *,'ERR=',MAXVAL(ABS(F-F0))
    END SUBROUTINE EX06
    
    SUBROUTINE PENRES1( N, A, B, C, D, E, F, X, R)
        INTEGER, INTENT(IN) :: N
        REAL(8), INTENT(IN) :: A(*), B(*), C(*), D(*), E(*), F(*), X(*)
        REAL(8), INTENT(OUT) :: R(*)
        INTEGER :: K
        DO K = 1, N
            IF (K == 1) THEN
                R(K) = C(K)*X(K) + D(K)*X(K+1) + E(K)*X(K+2) + A(1)*X(N-1) + B(K)*X(N) - F(K)
            ELSE IF (K == 2) THEN
                R(K) = B(K)*X(K-1) + C(K)*X(K) + D(K)*X(K+1) + E(K)*X(K+2) + A(K)*X(N) - F(K)            
            ELSE IF (K == N-1) THEN
                R(K) = E(K)*X(1) + A(K)*X(K-2) + B(K)*X(K-1) + C(K)*X(K) + D(K)*X(K+1)  - F(K)            
            ELSE IF (K == N) THEN
                R(K) = D(K)*X(1) + E(K)*X(2) + A(K)*X(K-2) + B(K)*X(K-1) + C(K)*X(K)   - F(K)            
            ELSE
                R(K) = A(K)*X(K-2) + B(K)*X(K-1) + C(K)*X(K) + D(K)*X(K+1) + E(K)*X(K+2) - F(K)                        
            END IF
        ENDDO
    END SUBROUTINE PENRES1    
!
    SUBROUTINE PENRES2( M, N, A, B, C, D, E, F, X, R)
        INTEGER, INTENT(IN) :: N
        REAL(8), INTENT(IN) :: A(M,M,*), B(M,M,*), C(M,M,*), D(M,M,*), E(M,M,*), F(M,*), X(M,*)
        REAL(8), INTENT(OUT) :: R(M,*)
        INTEGER :: K
        DO K = 1, N
            IF (K == 1) THEN
                R(:,K) = MATMUL(C(:,:,K),X(:,K)) + MATMUL(D(:,:,K),X(:,K+1)) &
                &   + MATMUL(E(:,:,K),X(:,K+2)) + MATMUL(A(:,:,K),X(:,N-1))  &
                &   + MATMUL(B(:,:,K),X(:,N)) - F(:,K)
            ELSE IF (K == 2) THEN
                R(:,K) = MATMUL(B(:,:,K),X(:,K-1)) + MATMUL(C(:,:,K),X(:,K)) &
                &   + MATMUL(D(:,:,K),X(:,K+1)) + MATMUL(E(:,:,K),X(:,K+2))  &
                &   + MATMUL(A(:,:,K),X(:,N)) - F(:,K)
            ELSE IF (K == N-1) THEN
                R(:,K) = MATMUL(A(:,:,K),X(:,K-2)) + MATMUL(B(:,:,K),X(:,K-1))&
                &   + MATMUL(C(:,:,K),X(:,K)) + MATMUL(D(:,:,K),X(:,K+1)) &
                &   + MATMUL(E(:,:,K),X(:,1))- F(:,K)            
            ELSE IF (K == N) THEN
                R(:,K) = MATMUL(A(:,:,K),X(:,K-2)) + MATMUL(B(:,:,K),X(:,K-1))&
                &   + MATMUL(C(:,:,K),X(:,K)) + MATMUL(D(:,:,K),X(:,1)) &
                &   + MATMUL(E(:,:,K),X(:,2))- F(:,K)            
            ELSE
                R(:,K) = MATMUL(A(:,:,K),X(:,K-2)) + MATMUL(B(:,:,K),X(:,K-1))&
                &   + MATMUL(C(:,:,K),X(:,K)) + MATMUL(D(:,:,K),X(:,K+1)) &
                &   + MATMUL(E(:,:,K),X(:,K+2))- F(:,K)
            END IF
        ENDDO
    END SUBROUTINE PENRES2       
!                
END PROGRAM    
