SUBROUTINE STLIO4(ELAM,LEFT,RIGHT,MATCHP,SIZE,TOL,BOUND_ROUT,PROB_DEF,PROB_SOL,LSOL)

! This routine computes the eigenfunctions of 4th order sturm liouville 
! problems, using continous orthonormalization.
!
! PARAMETERS :
! -----------
!
! *- ELAM : Complex(dble precision) 
!           On entry, ELAM is the eigenvalue that the user wants to use
!           for the computation. ELAM is unchanged on exit.
! *- LEFT : Real(dble precision)
!           On entry, LEFT is the left interval end of integration.
!           LEFT is unchanged on exit.
! *- RIGHT : Real(dble precision)
!           On entry, RIGHT is the right interval end of integration.
!           RIGHT is unchanged on exit.
! *- MATCHP : Real(dble precision)
!           On entry, MATCHP is the matching point the user wants to
!           define for the integration problem. 
!           MATCHP is unchanged on exit.
!           Constraint : LEFT < MATCHP < RIGHT
! *- SIZE : Integer
!           On entry, SIZE is approximately double of the size of solution
!           array that the user wants to obtain. This size will be adjust 
!           by the routine during the computation.
!           If the value supplied by the user is not big enough, a failure
!           message permit to warn the user, who can try with a bigger value.
!           SIZE is unchanged on exit.
! *- TOL : Real(dble precision)
!           On entry, TOL is the scalar tolerance defined by the user
!           for the integration. TOL is unchanged on exit.
!           Constraint : TOL > 10.0 * machine precision
! *- BOUND_ROUT : Subroutine, supplied by the user.
!           BOUND_ROUT gives the values of y and of its first derivative
!           at interval ends.
!           Its specification is :
!-----------------------------------------------------------------------
!           BOUND_ROUT(IEND,X,ELAM,U,V)
!
!      *- IEND : integer (= 0 or 1)
!             On entry, IEND indicates if the BC are needed at the left
!             interval end (IEND=0), or at the right interval end (IEND=1)
!             IEND is unchanged on exit.
!      *- X : Real (Dble precision)
!      *- ELAM : Complex (Dble precision)
!             The eigenvalue of the problem, which must be the same as
!             the value given to STLIO4.
!             ELAM is unchanged on exit.
!      *- U : Complex array (2,2)
!             On exit, U defines the upper part of the 4 by 2 matrix which
!             defines the BC.
!      *- V : Complex array (2,2)
!             On exit, V defineS the lower part of the 4 by 2 matrix which
!             defines the BC. 
!------------------------------------------------------------------------
! *- PROB_DEF :  Subroutine, supplied by the user.         
!           PROB_DEF gives the values of the differential equation 
!           coefficients P, Q3, Q2, Q1, Q0 which define the problem
!           for one value of X between LEFT and RIGHT
!           Its specification is :
!------------------------------------------------------------------------
!           PROB_DEF(X,ELAM,P,Q3,Q2,Q1,Q0)
!
!      *- X : Real (Dble precision)
!             On entry, X is the value at which the coeffients of the 
!             differential equation are needed.
!      *- ELAM : Complex (Dble precision)
!             The eigenvalue of the problem, which must be the same as
!             the value given to STLIO4.
!             ELAM is unchanged on exit.
!      *- P : Real (Dble precision)
!      *- Q3,Q2,Q1,Q0 : Complex (Dble precision)
!             On exit, P, Q3, Q2, Q1, Q0 are the differntial equation 
!             coefficients which defines the problem.
!------------------------------------------------------------------------
! *- PROB_SOL : Real array (0:2,1:SIZE)
!           On exit, PROB_SOL contains the computation results with in its first   
!           row the different times at which a solution has been computed, in its
!           second row the real part of the solution at this different times, and 
!           the last row the imaginary part of the solution at this times.
! *- LSOL : Integer
!           On exit, LSOL is the useful size of the solution array.
!           LSOL is approximately equal to SIZE/2.


IMPLICIT NONE

!....... Intent arguments .........!
Integer, intent(in) :: SIZE
Real(kind=kind(1.d0)), dimension(0:2,SIZE), intent(out) :: PROB_SOL
Real(kind=kind(1.d0)), intent(in) :: TOL, LEFT, RIGHT, MATCHP
Complex(kind=kind(1.d0)), intent(in) :: ELAM
Integer, intent(out) :: LSOL
EXTERNAL BOUND_ROUT,PROB_DEF

!....... Integer arguments ...........!
Integer :: LL, LR, SIZEP

!........ Real arguments & arrays .............!
Real(kind=kind(1.d0)), dimension(:), allocatable :: TIMER, TIMEL
Real(kind=kind(1.d0)) :: eps

!......... Complex arguments & arrays .......!
Complex(kind=kind(1.d0)), dimension(:,:), allocatable :: YSOLL, YSOLR
Complex(kind=kind(1.d0)), dimension(2) :: VECTL,VECTR
Complex(kind=kind(1.d0)), dimension(2,2) :: WRC,WLC

!.....Interface......!

INTERFACE

SUBROUTINE CALCEIGF(PI,PF,TOL,LEFT,ELAM,k,TIME,YSOL,SIZE,BOUND_ROUT,PROB_DEF)
  Integer, intent(in) :: SIZE 
  Real(kind=kind(1.d0)), intent(in) :: PI, PF, LEFT, TOL
  Complex(kind=kind(1.d0)), intent(in) :: ELAM
  Integer, intent(out) :: k
  Complex(kind=kind(1.d0)), dimension(1:12,1:SIZE), intent(out) :: YSOL 
  Real(kind=kind(1.d0)), dimension(1:SIZE), intent(out) :: TIME 
  Complex(kind=kind(1.d0)), dimension(1:2,1:2) :: U, V, W
  Complex(kind=kind(1.d0)), dimension(1:4,1:2) :: UV, T
  Complex(kind=kind(1.d0)), dimension(:), allocatable :: Y_USE, YNOW_USE
  Real(kind=kind(1.d0)) :: TNOW, eps, TCRIT,TOL, HMAX, TOUT, TINC
  Real(kind=kind(1.d0)),dimension(:), allocatable :: RWORK, ATOL, RTOL
  Real(kind=kind(1.d0)),dimension(:), allocatable :: YNOW
  Real(kind=kind(1.d0)) :: var, ceiling
  Integer :: NEQ, NEQG, HALFNEQ, LRWORK, IFAIL, IEND, LATOL, LRTOL, LIWORK, MAXSTP
  Integer, dimension(:), allocatable :: IWORK
  Integer :: IU, IV
  Integer :: i,j,z
  Character(1) :: STATEF
  Logical :: VECTOL, ONESTP, CRIT, ALTERG, SOPHST
  EXTERNAL PROB_DEF, BOUND_ROUT
END SUBROUTINE 

SUBROUTINE F(X,ELAM,PROB_DEF,Y,YP)
  Real(kind=kind(1.d0)), intent(in) :: X
  Real(kind=kind(1.d0)), dimension(24), intent(in) :: Y
  Real(kind=kind(1.d0)), dimension(24), intent(out) :: YP
  Complex(kind=kind(1.d0)), dimension(4,4) :: A
  Complex(kind=kind(1.d0)), dimension(12) :: YP_USE
  Complex(kind=kind(1.d0)), dimension(12) :: Y_USE
  Complex(kind=kind(1.d0)), intent(in) :: ELAM
  Complex(kind=kind(1.d0)), dimension(4,2) :: T, RES_T
  Complex(kind=kind(1.d0)), dimension(2,2) :: W, RES_W
  EXTERNAL PROB_DEF
END SUBROUTINE

SUBROUTINE ADJUST(YSOLL,YSOLR,LL,LR,SIZE,VECTL,VECTR,WRC,WLC,PROB_SOL,LSOL,TIMEL,TIMER)
  Integer, intent(in) :: SIZE, LL, LR
  Complex(kind=kind(1.d0)), dimension(1:12,1:SIZE), intent(in) :: YSOLL, YSOLR 
  Complex(kind=kind(1.d0)), dimension(2), intent(in) :: VECTL, VECTR
  Complex(kind=kind(1.d0)), dimension(2,2), intent(in) :: WRC, WLC 
  Real(kind=kind(1.d0)), dimension(SIZE), intent(in) :: TIMEL, TIMER
  Real(kind=kind(1.d0)), dimension(0:2,2*SIZE), intent(out) :: PROB_SOL
  Integer, intent(out) :: LSOL 
  Complex(kind=kind(1.d0)), dimension(12) :: Y_USE
  Complex(kind=kind(1.d0)), dimension(4,2) :: T_USE
  Complex(kind=kind(1.d0)), dimension(2,2) :: W_USE
  Complex(kind=kind(1.d0)), dimension(2) :: WC_VECT, W_WC_VECT
  Complex(kind=kind(1.d0)), dimension(4) :: SOL
  Complex(kind=kind(1.d0)) :: alpha,beta
  Character(1) :: NORM
  Integer :: i,j
END SUBROUTINE 

SUBROUTINE VECTFIND(YSOLL,YSOLR,LL,LR,SIZE,VECTL,VECTR,WLC,WRC)
  Integer, intent(in) :: SIZE, LL, LR
  Complex(kind=kind(1.d0)), dimension(1:12,1:SIZE), intent(in) :: YSOLL,YSOLR 
  Complex(kind=kind(1.d0)), dimension(2), intent(out) :: VECTL,VECTR
  Complex(kind=kind(1.d0)), dimension(2,2), intent(out) :: WRC, WLC 
  Complex(kind=kind(1.d0)), dimension(1:12) :: YL_USE,YR_USE
  Complex(kind=kind(1.d0)), dimension(4,2) :: TL,TR
  Complex(kind=kind(1.d0)), dimension(2,2) :: WR, WL, W_TRANS
  Complex(kind=kind(1.d0)), dimension(4,4) :: MATT 
  Complex(kind=kind(1.d0)), dimension(4) :: VAL_CONT
  Complex(kind=kind(1.d0)), dimension(4,4) :: VECT_CONT
  Complex(kind=kind(1.d0)), dimension(8) :: WORK
  Complex(kind=kind(1.d0)), dimension(2) :: WORKINV
  Complex(kind=kind(1.d0)) :: VAL
  Real(kind=kind(1.d0)), dimension(8) :: RWORK
  Real(kind=kind(1.d0)) :: norme1, norme2
  Integer :: i,j, INDEX
  Integer :: IFAIL, LWORK, LDV, LDA
  Integer :: LDW,  INFO, LWORKINV
  Integer(2) :: IPIV
END SUBROUTINE

SUBROUTINE ROUT_INTEG(F,ELAM,PROB_DEF,NEQ,TNOW,TOUT,YNOW,RWORK,LRWORK,IWORK,LIWORK,IFAIL)
  Integer, intent(in) :: NEQ, LRWORK, LIWORK
  Real(kind=kind(1.d0)), dimension(LRWORK), intent(inout) :: RWORK
  Integer, dimension(LIWORK), intent(inout) :: IWORK
  Integer, intent(out) :: IFAIL
  Real(kind=kind(1.d0)), intent(inout) :: TNOW, TOUT
  Real(kind=kind(1.d0)), dimension(NEQ), intent(inout) :: YNOW
  Complex(kind=kind(1.d0)), intent(in) :: ELAM
  EXTERNAL PROB_DEF, F
  Integer :: NEQG, IREVCM, YRVCM, YPRVCM, KGRVCM, i
  Logical :: ROOT
  Real(kind=kind(1.d0)) :: TRVCM, GRVCM
  Real(kind=kind(1.d0)), dimension(NEQ) :: Y, YP
END SUBROUTINE

SUBROUTINE FACTQR(A,m,n,Q,R)
  Complex(kind=kind(1.d0)), Dimension(1:m,1:n), intent(in) :: A 
  Complex(kind=kind(1.d0)), Dimension(1:m,1:n), intent(out) :: Q
  Complex(kind=kind(1.d0)), Dimension(1:n,1:n), intent(out) :: R
  Integer, intent(in) :: m, n
  Complex(kind=kind(1.d0)), Dimension(1:n,1:n) :: RTRANS 
  Complex(kind=kind(1.d0)), Dimension(:), allocatable :: TAU 
  Complex(kind=kind(1.d0)), Dimension(:), allocatable :: WORK
  Real(kind=kind(1.d0)), Dimension(:), allocatable :: RWORK
  Integer :: INFO, i, j, LDA, LWORK
  Integer, Dimension(:), allocatable :: JPVT
END SUBROUTINE

SUBROUTINE VECT_MAT(T,mt,nt,W,mw,nw,Z,nz)
  Integer, intent(in) :: mt, nt
  Integer, intent(in) :: mw, nw
  Integer, intent(in) :: nz
  Complex(kind=kind(1.d0)), Dimension(1:mt,1:nt), intent(out) :: T  
  Complex(kind=kind(1.d0)), Dimension(1:mw,1:nw), intent(out) :: W
  Complex(kind=kind(1.d0)), Dimension(1:nz), intent(in) :: Z
  Integer :: i,j
END SUBROUTINE

SUBROUTINE MAT_VECT(T,mt,nt,W,mw,nw,Z)
  Integer, intent(in) :: mt, nt
  Integer, intent(in) :: mw, nw
  Complex(kind=kind(1.d0)), Dimension(1:mt,1:nt), intent(in) :: T  
  Complex(kind=kind(1.d0)), Dimension(1:mw,1:nw), intent(in) :: W
  Complex(kind=kind(1.d0)), Dimension(1:24), intent(out) :: Z
  Integer :: i,j,k,l
END SUBROUTINE

SUBROUTINE REAL_COMP(AR,n,ACOMP)
  Integer, intent(in) :: n
  Complex(kind=kind(1.d0)), Dimension(1:n), intent(out) :: ACOMP
  Real(kind=kind(1.d0)), Dimension(1:2*n), intent(in):: AR 
  Integer :: i
END SUBROUTINE

SUBROUTINE COMP_REAL(ACOMP,n,AR)
  Integer, intent(in) :: n
  Complex(kind=kind(1.d0)), Dimension(1:n), intent(in) :: ACOMP
  Real(kind=kind(1.d0)), Dimension(1:2*n), intent(out):: AR 
  Integer :: i
END SUBROUTINE

SUBROUTINE CALC_U(T,W,A,RES_T,RES_W)
  Complex(kind=kind(1.d0)), Dimension(1:4,1:4), intent(in) :: A
  Complex(kind=kind(1.d0)), Dimension(1:4,1:2), intent(in) :: T
  Complex(kind=kind(1.d0)), Dimension(1:2,1:2), intent(in) :: W
  Complex(kind=kind(1.d0)), Dimension(1:2,1:2), intent(out) :: RES_W
  Complex(kind=kind(1.d0)), Dimension(1:4,1:2), intent(out) :: RES_T
  Complex(kind=kind(1.d0)), Dimension(1:2,1:2) :: U
  Complex(kind=kind(1.d0)), Dimension(1:4,1:2) :: A_T
  Complex(kind=kind(1.d0)), Dimension(1:2,1:2) :: TSTAR_A_T
  Complex(kind=kind(1.d0)), Dimension(1:2,1:2) :: TSTAR_T
  Complex(kind=kind(1.d0)), Dimension(1:2,1:2) :: TSTAR_T_INV
  Complex(kind=kind(1.d0)), Dimension(1:4,1:2) :: T_U
  Integer :: i,j
  Character(1) :: NORM, STAR
  Complex(kind=kind(1.d0)) :: alpha, beta
  Integer, dimension(2) :: IPIV
  Integer :: INFO, LWORK
  Complex(kind=kind(1.d0)), Dimension(:), allocatable :: WORK
  Complex(kind=kind(1.d0)), Dimension(1:2,1:2) :: ATRANS
END SUBROUTINE

SUBROUTINE CALC_A(X,ELAM,PROB_DEF,A)
  Complex(kind=kind(1.d0)), Dimension(1:4,1:4) :: A
  Complex(kind=kind(1.d0)),intent(in) :: ELAM
  Real(kind=kind(1.d0)),intent(in) :: X
  EXTERNAL PROB_DEF
  Complex(kind=kind(1.d0)) :: P,Q3,Q2,Q1,Q0
END SUBROUTINE

END INTERFACE

!.... PROGRAM ........ PROGRAM ........ PROGRAM ........ PROGRAM ........ PROGRAM .....!

!........ Allocate arrays ..........!

Allocate(YSOLL(1:12,1:SIZE))
Allocate(YSOLR(1:12,1:SIZE))
Allocate(TIMEL(1:SIZE))
Allocate(TIMER(1:SIZE))

!....... Executable statements ...........!

SIZEP=SIZE/2

Call CALCEIGF(LEFT,MATCHP,TOL,LEFT,ELAM,LL,TIMEL,YSOLL,SIZEP,BOUND_ROUT,PROB_DEF)
Call CALCEIGF(RIGHT,MATCHP,TOL,LEFT,ELAM,LR,TIMER,YSOLR,SIZEP,BOUND_ROUT,PROB_DEF)
Call VECTFIND(YSOLL,YSOLR,LL,LR,SIZEP,VECTL,VECTR,WLC,WRC)
Call ADJUST(YSOLL,YSOLR,LL,LR,SIZEP,VECTL,VECTR,WRC,WLC,PROB_SOL,LSOL,TIMEL,TIMER)

END SUBROUTINE STLIO4

!---------------------------------------------------------------------------------

SUBROUTINE ADJUST(YSOLL,YSOLR,LL,LR,SIZE,VECTL,VECTR,WRC,WLC,PROB_SOL,LSOL,TIMEL,TIMER)

! This routine is the final routine which adjusts the two parts of the solution

IMPLICIT NONE

!.......... Intent arguments ........!
Integer, intent(in) :: SIZE, LL, LR
Complex(kind=kind(1.d0)), dimension(1:12,1:SIZE), intent(in) :: YSOLL, YSOLR 
Complex(kind=kind(1.d0)), dimension(2), intent(in) :: VECTL, VECTR
Complex(kind=kind(1.d0)), dimension(2,2), intent(in) :: WRC, WLC 
Real(kind=kind(1.d0)), dimension(SIZE), intent(in) :: TIMEL, TIMER
Real(kind=kind(1.d0)), dimension(0:2,2*SIZE), intent(out) :: PROB_SOL
Integer, intent(out) :: LSOL 

!........ Local arguments ...........!
Complex(kind=kind(1.d0)), dimension(12) :: Y_USE
Complex(kind=kind(1.d0)), dimension(4,2) :: T_USE
Complex(kind=kind(1.d0)), dimension(2,2) :: W_USE
Complex(kind=kind(1.d0)), dimension(2) :: WC_VECT, W_WC_VECT
Complex(kind=kind(1.d0)), dimension(4) :: SOL
Complex(kind=kind(1.d0)) :: alpha,beta
Character(1) :: NORM
Integer :: i,j

!...... Executable statements ......!

NORM = 'N'
alpha = DCMPLX(1.0,0.0)
beta = DCMPLX(0.0,0.0)

LSOL = LL+LR

DO i=1,LL
   DO j=1,12
      Y_USE(j) = YSOLL(j,i)
   END DO
   Call VECT_MAT(T_USE,4,2,W_USE,2,2,Y_USE,12)
   Call F06ZAF(NORM,NORM,2,1,2,alpha,WLC,2,VECTL,2,beta,WC_VECT,2)
   Call F06ZAF(NORM,NORM,2,1,2,alpha,W_USE,2,WC_VECT,2,beta,W_WC_VECT,2)
   Call F06ZAF(NORM,NORM,4,1,2,alpha,T_USE,4,W_WC_VECT,2,beta,SOL,4)
   PROB_SOL(1,i) = Dble(SOL(1))
   PROB_SOL(2,i) = Dimag(SOL(1))
   PROB_SOL(0,i) = TIMEL(i) 
END DO

DO i=1,LR
   DO j=1,12
      Y_USE(j) = YSOLR(j,i)
   END DO
   Call VECT_MAT(T_USE,4,2,W_USE,2,2,Y_USE,12)
   Call F06ZAF(NORM,NORM,2,1,2,alpha,WRC,2,VECTR,2,beta,WC_VECT,2)
   Call F06ZAF(NORM,NORM,2,1,2,alpha,W_USE,2,WC_VECT,2,beta,W_WC_VECT,2)
   Call F06ZAF(NORM,NORM,4,1,2,alpha,T_USE,4,W_WC_VECT,2,beta,SOL,4) 
   PROB_SOL(1,LL+i) = Dble(SOL(1))
   PROB_SOL(2,LL+i) = Dimag(SOL(1))
   PROB_SOL(0,LL+i) = TIMER(i)
END DO

END SUBROUTINE

!---------------------------------------------------------------------------------

SUBROUTINE VECTFIND(YSOLL,YSOLR,LL,LR,SIZE,VECTL,VECTR,WLC,WRC)

! This routine is used in order to find how the adjust the two parts
! of the solution, searching for the eigenvector corresponding to the 
! the 0 eigenvalue of a particular matrix.

IMPLICIT NONE

!.......... Intent arguments ........!
Integer, intent(in) :: SIZE, LL, LR
Complex(kind=kind(1.d0)), dimension(1:12,1:SIZE), intent(in) :: YSOLL,YSOLR 
Complex(kind=kind(1.d0)), dimension(2), intent(out) :: VECTL,VECTR
Complex(kind=kind(1.d0)), dimension(2,2), intent(out) :: WRC, WLC 

!......... Local arguments ........!
Complex(kind=kind(1.d0)), dimension(1:12) :: YL_USE,YR_USE
Complex(kind=kind(1.d0)), dimension(4,2) :: TL,TR
Complex(kind=kind(1.d0)), dimension(2,2) :: WR, WL, W_TRANS
Complex(kind=kind(1.d0)), dimension(4,4) :: MATT 
Complex(kind=kind(1.d0)), dimension(4) :: VAL_CONT
Complex(kind=kind(1.d0)), dimension(4,4) :: VECT_CONT
Complex(kind=kind(1.d0)), dimension(8) :: WORK
Complex(kind=kind(1.d0)), dimension(2) :: WORKINV
Complex(kind=kind(1.d0)) :: VAL

Real(kind=kind(1.d0)), dimension(8) :: RWORK
Real(kind=kind(1.d0)) :: norme1, norme2

Integer :: i,j, INDEX
Integer :: IFAIL, LWORK, LDV, LDA
Integer :: LDW,  INFO, LWORKINV
Integer(2) :: IPIV

!....... Executable statements .....!

DO i=1,12 
   YL_USE(i) = YSOLL(i,LL)
   YR_USE(i) = YSOLR(i,1)
END DO

Call VECT_MAT(TR,4,2,WR,2,2,YR_USE,12)
Call VECT_MAT(TL,4,2,WL,2,2,YL_USE,12)

LDW = 2
INFO = 0
W_TRANS = WR
LWORKINV = 2
Call F07ARF(LDW,LDW,W_TRANS,LDW,IPIV,INFO)
IF (INFO /= 0) Write(*,*) 'Problem1 with the matrix inversion, INFO =',INFO
Call F07AWF(LDW,W_TRANS,LDW,IPIV,WORKINV,LWORKINV,INFO)
IF (INFO /= 0) Write(*,*) 'Problem2 with the matrix inversion, INFO =',INFO
WRC = W_TRANS

W_TRANS = WL
IPIV = 0
WORKINV = DCMPLX(0.0,0.0)
Call F07ARF(LDW,LDW,W_TRANS,LDW,IPIV,INFO)
IF (INFO /= 0) Write(*,*) 'Problem1 with the matrix inversion, INFO =',INFO
Call F07AWF(LDW,W_TRANS,LDW,IPIV,WORKINV,LWORKINV,INFO)
IF (INFO /= 0) Write(*,*) 'Problem2 with the matrix inversion, INFO =',INFO
WLC = W_TRANS

DO i=1,4
   DO j=1,2
      MATT(i,j) = TL(i,j)
      MATT(i,j+2) = -TR(i,j)      
   END DO
END DO

LDA = 4
LDV = 4
LWORK = 2*LDA
IFAIL = 0
Call F02GBF('V',LDA,MATT,LDA,VAL_CONT,VECT_CONT,LDV,RWORK,WORK,LWORK,IFAIL)

IF (IFAIL /= 0)  Write(*,*) 'Problem with the eigenvector computation'

VAL = VAL_CONT(1)
INDEX = 1
norme1 = (Dble(VAL)**2+Dimag(VAL)**2)**(0.5)
DO i=2,4
   norme2 = (Dble(VAL_CONT(i))**2+Dimag(VAL_CONT(i))**2)**(0.5)
   IF ( norme2 <= norme1 ) THEN
      VAL = VAL_CONT(i)
      INDEX = i
      norme1 = norme2
   END IF
END DO

DO i=1,2
   VECTL(i) = VECT_CONT(i,INDEX)
   VECTR(i) = VECT_CONT(i+2,INDEX)
END DO

END SUBROUTINE

!----------------------------------------------------------------------------------

SUBROUTINE CALCEIGF(PI,PF,TOL,LEFT,ELAM,k,TIME,YSOL,SIZE,BOUND_ROUT,PROB_DEF)

! This routine computes the eigenfunction F
! NOTE : Y is the solution of the fourth order problem :
!                            L(lambda)Y = 0
!
! PARAMETERS :
! ------------
! *- PI : real (input) : on entry, PI is the point where starts 
!                        the integration. (unchanged on exit)
! *- PF : real (output) : on entry, PF is the point where ends 
!                         the integration. (unchanged on exit)
! *- TOL : real (input) : on entry, TOL is the error tolerance we have 
!                         to define for the NAG routines.
!                         (unchanged on exit)
! *- LEFT : real (input) : on entry, LEFT is a reference point, which permit to define
!                       the problem (boundary conditions, the sense of the 
!                       integration). (unchanged on exit)
! *- ELAM  : real (input) : on entry, Lambda is the parameter needed to define
!                            the operator L. The different possible values of
!                            LAMBDA are the eigenvalues given by the routine
!                            EIGVAL. (unchanged on exit)
! *- k : integer (output) : k give us a very important information, the quantity 
!                           of datas stored in the different arrays. So, the 
!                           parameter SIZE has to be greater than k. Otherwise, 
!                           the routine fails with an error message to warn the 
!                           user.
! *- TIME : real array(SIZE) , output : TIME table contains the time values at which
!                            a solution has been computed
! *- YSOL : complex array(12,SIZE), output : YSOL contains all the information about
!                            The computation of the matrix T and W such as :
!                                      Y = T * W
!                            This informations are in the vector YSOL(:,i) for the ith
!                            step of the integration
! *- SIZE : integer (input) : SIZE is the YSOL and TIME size chosen by the user before 
!                             to call this routine.
!                             WARNING : it is not the number of datas in output
! *- BOUND_ROUT : external routine (input) : BOUND_ROUT is the routine supplied 
!                    by the user which defines the boudary conditions at left
!                    and right interval ends.
! *- PROB_DEF : external routine (input) : PROB_DEF is the routine supplied 
!                    by the user which defines the differential equation 
!                    coefficients.

IMPLICIT NONE

!..... Common arguments ........!
Integer, intent(in) :: SIZE 

!.......... Intent arguments ........!
Real(kind=kind(1.d0)), intent(in) :: PI, PF, LEFT, TOL
Complex(kind=kind(1.d0)), intent(in) :: ELAM
Integer, intent(out) :: k
Complex(kind=kind(1.d0)), dimension(1:12,1:SIZE), intent(out) :: YSOL 
Real(kind=kind(1.d0)), dimension(1:SIZE), intent(out) :: TIME 

!..... Local complex arrays & arguments .....!
Complex(kind=kind(1.d0)), dimension(1:2,1:2) :: U, V, W
Complex(kind=kind(1.d0)), dimension(1:4,1:2) :: UV, T
Complex(kind=kind(1.d0)), dimension(:), allocatable :: Y_USE, YNOW_USE

!..... Local real arrays & arguments ........!
Real(kind=kind(1.d0)) :: TNOW, eps, TCRIT,TOL, HMAX, TOUT, TINC
Real(kind=kind(1.d0)),dimension(:), allocatable :: RWORK, ATOL, RTOL
Real(kind=kind(1.d0)),dimension(:), allocatable :: YNOW
Real(kind=kind(1.d0)) :: var, ceiling

!........ Integer arguments .........!
Integer :: NEQ, NEQG, HALFNEQ, LRWORK, IFAIL, IEND, LATOL, LRTOL, LIWORK, MAXSTP
Integer, dimension(:), allocatable :: IWORK
Integer :: IU, IV
Integer :: i,j,z
 
!......... Other arguments ..........!
Character(1) :: STATEF
Logical :: VECTOL, ONESTP, CRIT, ALTERG, SOPHST
EXTERNAL PROB_DEF, BOUND_ROUT, F

!....... Executable statements .........!

eps = TOL
NEQ = 24
HALFNEQ = NEQ/2
LRWORK = 32*NEQ
IFAIL = -1
ceiling = 10.0
STATEF = 'S'
LATOL = 1
LRTOL = 1
VECTOL = .false.
ONESTP = .false.
CRIT = .true.
TCRIT = PF
HMAX = 0.0
MAXSTP = 0
NEQG = 0
ALTERG = .false.
SOPHST = .false.
LRWORK = 23*(1+NEQ)+4
LIWORK = 21 + NEQG
TINC = 2.0 * (PF-PI) / SIZE

Allocate(YNOW(NEQ))
Allocate(Y_USE(HALFNEQ))
Allocate(YNOW_USE(HALFNEQ))
Allocate(ATOL(LATOL))
Allocate(RTOL(LRTOL))
Allocate(RWORK(LRWORK))
Allocate(IWORK(LIWORK))

ATOL(1) = TOL
RTOL(1) = TOL

IF ((PI-LEFT) <= eps) THEN 
   IEND = 0
ELSE
   IEND = 1
END IF

Call BOUND_ROUT(IEND,PI,ELAM,U,V)
DO i=1,2
   DO j=1,2
      UV(i,j) = U(i,j)
      UV(i+2,j) = V(i,j)
   END DO
END DO

Call FACTQR(UV,4,2,T,W)
Call MAT_VECT(T,4,2,W,2,2,Y_USE)
Call COMP_REAL(Y_USE,12,YNOW)

k=1

TIME(k) = PI

DO i=1,12
   YSOL(i,k) = Y_USE(i)
END DO

IFAIL = 0

Call D02QWF(STATEF,NEQ,VECTOL,ATOL,LATOL,RTOL,LRTOL,ONESTP,CRIT,TCRIT,HMAX,MAXSTP,NEQG,&
     &    ALTERG,SOPHST,RWORK,LRWORK,IWORK,LIWORK,IFAIL)

TNOW =  PI

DO WHILE (abs(TNOW-PF)>eps)
   
   IF (k >= SIZE) THEN
      Write(*,*) 'The table dimension is not big enough, try with a bigger value.'
      STOP
   ELSE
      IFAIL = -1
      TOUT = TNOW  + TINC
      IF ( PI < PF .and. TOUT > TCRIT) THEN
         TOUT = PF
      ELSE IF ( PI > PF .and. TOUT < TCRIT) THEN
         TOUT = PF
      END IF
      
      Call ROUT_INTEG(F,ELAM,PROB_DEF,NEQ,TNOW,TOUT,YNOW,RWORK,LRWORK,IWORK,LIWORK,IFAIL)
      
      IF (IFAIL==0) THEN
         
         IF ((PI-LEFT) <= eps) THEN 
            TIME(k+1) = TNOW
            Call REAL_COMP(YNOW,12,YNOW_USE)
            DO i=1,12
               YSOL(i,k+1) = YNOW_USE(i)
            END DO
         ELSE
            DO i=1,k
               j = k-i+1
               TIME(j+1) = TIME(j)
               DO z=1,12
                  YSOL(z,j+1) = YSOL(z,j)
               END DO
            END DO
            Call REAL_COMP(YNOW,12,YNOW_USE)               
            TIME(1) = TNOW
            DO i=1,12
               YSOL(i,1) = YNOW_USE(i)
            END DO
         END IF
         
         var = max(abs(YNOW(17)),abs(YNOW(18)),abs(YNOW(19)),abs(YNOW(20)))
         var = max(var,abs(YNOW(21)),abs(YNOW(22)),abs(YNOW(23)),abs(YNOW(24)))
         
         IF ( var >= ceiling ) THEN           
            DO i=1,k+1 
               DO j=9,12
                  YSOL(j,i) = YSOL(j,i)/ceiling
               END DO
            END DO
            DO i=1,16
               YNOW(i) = YNOW(i)
            END DO
            DO i=17,24
               YNOW(i) = YNOW(i)/ceiling
            END DO
            IF (abs(TNOW-PF)>=eps) THEN
               STATEF = 'R'
               Call D02QWF(STATEF,NEQ,VECTOL,ATOL,LATOL,RTOL,LRTOL,ONESTP,CRIT,TCRIT, &
                    & HMAX,MAXSTP,NEQG,ALTERG,SOPHST,RWORK,LRWORK,IWORK,LIWORK,IFAIL)
            END IF
         ELSE
            STATEF = 'C'
            Call D02QWF(STATEF,NEQ,VECTOL,ATOL,LATOL,RTOL,LRTOL,ONESTP,CRIT,TCRIT, &
                 & HMAX,MAXSTP,NEQG,ALTERG,SOPHST,RWORK,LRWORK,IWORK,LIWORK,IFAIL)           
         END IF
      ELSE
         Write(*,*)'IFAIL =', IFAIL
         STOP
      END IF
   END IF   
   k = k+1
END DO

END SUBROUTINE

!--------------------------------------------------------------------------------

SUBROUTINE ROUT_INTEG(F,ELAM,PROB_DEF,NEQ,TNOW,TOUT,YNOW,RWORK,LRWORK,IWORK,LIWORK,IFAIL)

! This routine is my own integrating routine which uses the D02QGF nag routine,
! with reverse communication.

IMPLICIT NONE

!...... Intent arguments ........!
Integer, intent(in) :: NEQ, LRWORK, LIWORK
Real(kind=kind(1.d0)), dimension(LRWORK), intent(inout) :: RWORK
Integer, dimension(LIWORK), intent(inout) :: IWORK
Integer, intent(out) :: IFAIL
Real(kind=kind(1.d0)), intent(inout) :: TNOW, TOUT
Real(kind=kind(1.d0)), dimension(NEQ), intent(inout) :: YNOW
Complex(kind=kind(1.d0)), intent(in) :: ELAM
EXTERNAL PROB_DEF, F

!....... Local arguments ..........!
Integer :: NEQG, IREVCM, YRVCM, YPRVCM, KGRVCM, i
Logical :: ROOT
Real(kind=kind(1.d0)) :: TRVCM, GRVCM
Real(kind=kind(1.d0)), dimension(NEQ) :: Y, YP

!........ Executable statements ........!

NEQG = 0
ROOT = .false.
IREVCM = 0
IFAIL = -1

Call D02QGF(NEQ,TNOW,YNOW,TOUT,NEQG,ROOT,IREVCM,TRVCM,YRVCM,YPRVCM,GRVCM,KGRVCM,RWORK,LRWORK,IWORK,LIWORK,IFAIL)

DO WHILE (IREVCM > 0 .and. IREVCM < 9)
   
   IF (IREVCM < 8) THEN
      IF (YRVCM == 0) THEN
         Y = YNOW
      ELSE
         DO i=1,NEQ
            Y(i) = RWORK(YRVCM+i-1)
         END DO
      END IF
      Call F(TRVCM,ELAM,PROB_DEF,Y,YP)
      DO i=1,NEQ
         RWORK(YPRVCM+i-1) = YP(i)
      END DO
   END IF
   
   IFAIL = -1
   Call D02QGF(NEQ,TNOW,YNOW,TOUT,NEQG,ROOT,IREVCM,TRVCM,YRVCM,YPRVCM,GRVCM,KGRVCM,RWORK,LRWORK,IWORK,LIWORK,IFAIL)
   
END DO

END SUBROUTINE 

!---------------------------------------------------------------------------------

SUBROUTINE F(X,ELAM,PROB_DEF,Y,YP)

! This routine is the routine needed to evaluate y prime

IMPLICIT NONE

!..... Intent arguments .......!
Real(kind=kind(1.d0)), intent(in) :: X
Real(kind=kind(1.d0)), dimension(24), intent(in) :: Y
Real(kind=kind(1.d0)), dimension(24), intent(out) :: YP

!...... Local arguments ........!
Complex(kind=kind(1.d0)), dimension(12) :: YP_USE
Complex(kind=kind(1.d0)), dimension(12) :: Y_USE

!.... Common arguments .........!
Complex(kind=kind(1.d0)), intent(in) :: ELAM
Complex(kind=kind(1.d0)), dimension(4,2) :: T, RES_T
Complex(kind=kind(1.d0)), dimension(2,2) :: W, RES_W
Complex(kind=kind(1.d0)), dimension(4,4) :: A
EXTERNAL PROB_DEF

!........ Executable statements  .........!

Call REAL_COMP(Y,12,Y_USE)
Call VECT_MAT(T,4,2,W,2,2,Y_USE,12)
Call CALC_A(X,ELAM,PROB_DEF,A)
Call CALC_U(T,W,A,RES_T,RES_W)
Call MAT_VECT(RES_T,4,2,RES_W,2,2,YP_USE)
Call COMP_REAL(YP_USE,12,YP)

END SUBROUTINE

!---------------------------------------------------------------------------------

SUBROUTINE FACTQR(A,m,n,Q,R)

! On entry : A(m,n) is the matrix you want to factorize 
!            m,n are the A dimensions
! On exit :  Q(m,n) and R(n,n) are the results of the factorization

IMPLICIT NONE

!......... Intent arguments ..........!
Complex(kind=kind(1.d0)), Dimension(1:m,1:n), intent(in) :: A 
Complex(kind=kind(1.d0)), Dimension(1:m,1:n), intent(out) :: Q
Complex(kind=kind(1.d0)), Dimension(1:n,1:n), intent(out) :: R
Integer, intent(in) :: m, n

!......... Local arguments ...........!
Complex(kind=kind(1.d0)), Dimension(1:n,1:n) :: RTRANS 
Complex(kind=kind(1.d0)), Dimension(:), allocatable :: TAU 
Complex(kind=kind(1.d0)), Dimension(:), allocatable :: WORK
Real(kind=kind(1.d0)), Dimension(:), allocatable :: RWORK
Integer :: INFO, i, j, LDA, LWORK
Integer, Dimension(:), allocatable :: JPVT

!......... Allocate arrays ...........!
Allocate(TAU(max(m,n)))
Allocate(WORK(3*n))
Allocate(RWORK(2*n))
Allocate(JPVT(n))
 
!....... Executable statements .........!

LDA = m
LWORK = 4*max(m,n)

Call F08ASF(m,n,A,LDA,TAU,WORK,LWORK,INFO)

DO i=1,n
   DO j=1,n
      IF (j >=i) THEN
         R(i,j) = A(i,j)
      ELSE
         R(i,j) = DCMPLX(0.0,0.0)
      END IF
   END DO
END DO

Call F08ATF(m,n,n,A,LDA,TAU,WORK,LWORK,INFO)

IF (INFO /= 0) Write(*,*) 'There is a QR factorisation problem' 

DO i=1,m
   DO j=1,n
      Q(i,j) = A(i,j)
   END DO
END DO

END SUBROUTINE

!-----------------------------------------------------------------------------------

SUBROUTINE MAT_VECT(T,mt,nt,W,mw,nw,Z)

! This routine transforms the two matrix T and W in one single vector Z in this way :
! 
!          (  Z(1)    Z(2)         .........           Z(nt)       )
!          (    .    .                                             )
!          (    .        .                                         )
! T(i,j) = (    .            Z((i-1)*mt+j*nt)                      )
!          (    .                                 .                )
!          (    .                                    .             )
!          (  Z((mt-1)*nt+1)         .........          Z(mt*nt)   )
!
!
!
!
!          (  Z(mt*nt+1)    Z(mt*nt+2)    .........     Z(mt*nt+nw)       )
!          (    .    .                                                    )
!          (    .        .                                                )
! T(i,j) = (    .            Z(mt*nt+(i-1)*mw+j*nw)                       )
!          (    .                                 .                       )
!          (    .                                     .                   )
!          (  Z(mt*nt+(mw-1)*nw+1)        .........      Z(mt*nt+mw*nw)   )
!
!
! PARAMETERS : T and W are the matrix you want to tranform into a vector
! with T(mt,nt), W(mw,nw) and Z(mt*nt+mw*nw)

IMPLICIT NONE

!......... Intent arguments ..........!
Integer, intent(in) :: mt, nt
Integer, intent(in) :: mw, nw
Complex(kind=kind(1.d0)), Dimension(1:mt,1:nt), intent(in) :: T  
Complex(kind=kind(1.d0)), Dimension(1:mw,1:nw), intent(in) :: W
Complex(kind=kind(1.d0)), Dimension(1:24), intent(out) :: Z

!........ Local arguments .............!
Integer :: i,j,k,l

!....... Executable statements .......!

DO k=1,(mt*nt+mw*nw)
   
   IF (k <= mt*nt) THEN
      i = (k-1)/nt + 1   
      j = k-(i-1)*nt
      Z(k) = T(i,j)
   ELSE
      l = k-mt*nt
      i = (l-1)/nw + 1
      j = l-(i-1)*nw
      Z(k) = W(i,j)
   END IF
   
END DO

END SUBROUTINE 

!-----------------------------------------------------------------------------------

SUBROUTINE VECT_MAT(T,mt,nt,W,mw,nw,Z,nz)

! This routine transforms the vector Z in two matrix T and W in this way :
! 
!          (  Z(1)    Z(2)         .........           Z(nt)       )
!          (    .    .                                             )
!          (    .        .                                         )
! T(i,j) = (    .                Z((i-1)*mt+j)                     )
!          (    .                                 .                )
!          (    .                                    .             )
!          (  Z((mt-1)*nt+1)         .........          Z(mt*nt)   )
!
!
!
!
!          (  Z(mt*nt+1)    Z(mt*nt+2)    .........     Z(mt*nt+nw)       )
!          (    .    .                                                    )
!          (    .        .                                                )
! T(i,j) = (    .             Z((i-1)*mw+j + mt*nt)                       )
!          (    .                                 .                       )
!          (    .                                     .                   )
!          (  Z(mt*nt+(mw-1)*nw+1)        .........      Z(mt*nt+mw*nw)   )
!
!
! PARAMETERS : T and W are the matrix you want to tranform into a vector
! with T(mt,nt), W(mw,nw) and Z(nz)

IMPLICIT NONE
  
!......... Intent arguments ..........!
Integer, intent(in) :: mt, nt
Integer, intent(in) :: mw, nw
Integer, intent(in) :: nz
Complex(kind=kind(1.d0)), Dimension(1:mt,1:nt), intent(out) :: T  
Complex(kind=kind(1.d0)), Dimension(1:mw,1:nw), intent(out) :: W
Complex(kind=kind(1.d0)), Dimension(1:nz), intent(in) :: Z

!........ Local arguments .............!
Integer :: i,j

!....... Executable statements .......!

IF ( nz /= mt*nt+mw*nw ) Write(*,*) 'Vector-matrix transformation impossible.'

DO i=1,mt
   DO j=1,nt
      T(i,j) = Z((i-1)*nt+j)
   END DO
END DO

DO i=1,mw
   DO j=1,nw
      W(i,j) = Z((i-1)*nw+j + mt*nt)
   END DO
END DO

END SUBROUTINE 

!-----------------------------------------------------------------------------------

SUBROUTINE REAL_COMP(AR,n,ACOMP)

! This routine converts the real vector AR(2n) into a complex vector ACOMP(n)
! with :
!      ACOMP(i) = (AR(2i-1),AR(2i))

IMPLICIT NONE

!.......... Intent arguments .............!
Integer, intent(in) :: n
Complex(kind=kind(1.d0)), Dimension(1:n), intent(out) :: ACOMP
Real(kind=kind(1.d0)), Dimension(1:2*n), intent(in):: AR 

!...... Local arguments ...........!
Integer :: i

!...... Executable statements .......!

DO i=1,n
   ACOMP(i) = DCMPLX(AR(2*i-1),AR(2*i))
END DO

END SUBROUTINE

!-----------------------------------------------------------------------------------

SUBROUTINE COMP_REAL(ACOMP,n,AR)

! This routine converts the complex vector ACOMP(n) into a real vector AR(2n)  
! with the real part of the first element, the imaginary part of the first
! element, the real part of the second element, etc

IMPLICIT NONE

!.......... Intent arguments .............!
Integer, intent(in) :: n
Complex(kind=kind(1.d0)), Dimension(1:n), intent(in) :: ACOMP
Real(kind=kind(1.d0)), Dimension(1:2*n), intent(out):: AR 

!...... Local arguments ...........!
Integer :: i

!...... Executable statements .......!

DO i=1,n
   AR(2*i-1) = DBLE(ACOMP(i))
   AR(2*i) = DIMAG(ACOMP(i))
END DO

END SUBROUTINE

!----------------------------------------------------------------------------

SUBROUTINE CALC_U(T,W,A,RES_T,RES_W)

! This routine calculates the U matrix which defines the two differential 
! equations :
!                    T' = A*T-TU
!                    W' = U*W
! With Y = T*W

IMPLICIT NONE

!....... Intent arguments ........!
Complex(kind=kind(1.d0)), Dimension(1:4,1:4), intent(in) :: A
Complex(kind=kind(1.d0)), Dimension(1:4,1:2), intent(in) :: T
Complex(kind=kind(1.d0)), Dimension(1:2,1:2), intent(in) :: W
Complex(kind=kind(1.d0)), Dimension(1:2,1:2), intent(out) :: RES_W
Complex(kind=kind(1.d0)), Dimension(1:4,1:2), intent(out) :: RES_T

!........ Local arguments ........!
Complex(kind=kind(1.d0)), Dimension(1:2,1:2) :: U
Complex(kind=kind(1.d0)), Dimension(1:4,1:2) :: A_T
Complex(kind=kind(1.d0)), Dimension(1:2,1:2) :: TSTAR_A_T
Complex(kind=kind(1.d0)), Dimension(1:2,1:2) :: TSTAR_T
Complex(kind=kind(1.d0)), Dimension(1:2,1:2) :: TSTAR_T_INV
Complex(kind=kind(1.d0)), Dimension(1:4,1:2) :: T_U
Integer :: i,j
Character(1) :: NORM, STAR
Complex(kind=kind(1.d0)) :: alpha, beta

!... Arguments for the inversion ........!
Integer, dimension(2) :: IPIV
Integer :: INFO, LWORK
Complex(kind=kind(1.d0)), Dimension(:), allocatable :: WORK
Complex(kind=kind(1.d0)), Dimension(1:2,1:2) :: ATRANS

!...... Executable statements ......!

NORM = 'N'
STAR = 'C'
alpha = DCMPLX(1.0,0.0)
beta = DCMPLX(0.0,0.0)
LWORK = 16
Allocate(WORK(LWORK))

Call F06ZAF(NORM,NORM,4,2,4,alpha,A,4,T,4,beta,A_T,4)
Call F06ZAF(STAR,NORM,2,2,4,alpha,T,4,A_T,4,beta,TSTAR_A_T,2)
Call F06ZAF(STAR,NORM,2,2,4,alpha,T,4,T,4,beta,TSTAR_T,2)

ATRANS = TSTAR_T
Call F07ARF(2,2,ATRANS,2,IPIV,INFO)

IF (INFO /= 0) THEN
   Write(*,*) 'Problem1 with the matrix inversion, INFO =',INFO
   DO i=1,2
      Write(*,*) (TSTAR_T(i,j),j=1,2)
   END DO
   STOP
END IF

Call F07AWF(2,ATRANS,2,IPIV,WORK,LWORK,INFO)
IF (INFO /= 0) Write(*,*) 'Problem2 with the matrix inversion, INFO =',INFO
TSTAR_T_INV = ATRANS

Call F06ZAF(NORM,NORM,2,2,2,alpha,TSTAR_T_INV,2,TSTAR_A_T,2,beta,U,2)
Call F06ZAF(NORM,NORM,4,2,2,alpha,T,4,U,2,beta,T_U,4)
Call F06ZAF(NORM,NORM,2,2,2,alpha,U,2,W,2,beta,RES_W,2)

DO i=1,4
   DO j=1,2
      RES_T(i,j) = A_T(i,j) - T_U(i,j)
   END DO
END DO

END SUBROUTINE

!----------------------------------------------------------------------------

SUBROUTINE CALC_A(X,ELAM,PROB_DEF,A)

! This routine calculates the A matrix of the differial problem Y' = A * Y
! PARAMETERS :
!      - X is the point at which the routine calculates the matrix
!      - ELAM is the eigenvalue the user has chosen for the problem
!      - PROB_DEF is the routine supplied by the user to define the
!        differential equation coefficients
!      - A (output) is the result of the calculation

IMPLICIT NONE

!....... Intent arguments ........!
Complex(kind=kind(1.d0)), Dimension(1:4,1:4) :: A
Complex(kind=kind(1.d0)),intent(in) :: ELAM
Real(kind=kind(1.d0)),intent(in) :: X
EXTERNAL PROB_DEF

!....... Local arguments .........!
Complex(kind=kind(1.d0)) :: P,Q3,Q2,Q1,Q0

!....... Executable statements ....!

Call PROB_DEF(X,ELAM,P,Q3,Q2,Q1,Q0)

A=0

A(1,2) = DCMPLX(1.0,0.0)
A(2,4) = - 1.0 / P
A(3,1) = - Q0
A(3,2) = - Q1
A(4,2) = Q2 
A(4,3) = DCMPLX(-1.0,0.0)
A(4,4) = - Q3 / P

END SUBROUTINE

!----------------------------------------------------------------------------