!****************************************************
!*   Find all roots of a complex polynomial using   *
!*   Laguerre formulation in complex domain.        *
!* ------------------------------------------------ *
!* SAMPLE RUN:                                      *
!* (Find roots of complex polynomial:               *
!*  (2.5 + I) X3 + (7.5 - 12 I) X2 + (-3.75 + 0.4 I *
!*  ) X + (4.25 - I)  )                             *
!*                                                  *
!*  Error code =         0                          *
!*                                                  *
!*  Complex roots are:                              *
!*                                                  *
!* (-1.079156,4.904068)                             *
!* (-0.1118458,0.6680070)                           *
!* (0.2599674,-0.3996614)                           *
!*                                                  *
!* ------------------------------------------------ *
!* Ref.: From Numath Library By Tuan Dang Trong in  *
!*       Fortran 77 [BIBLI 18].                     *
!*                                                  *
!*                F90 Release By J-P Moreau, Paris. *
!*                        (www.jpmoreau.fr)         *
!****************************************************   
PROGRAM TEST_LAGUERRE
parameter(N=3)  

COMPLEX  A(0:N), RAC(0:N), WORK(0:3*N)
integer IMP, ITMAX
real*8 EPS,EPS2

! Example #1
  A(0)= CMPLX(2.5d0,1.d0)
  A(1)= CMPLX(7.5d0,-12.d0)
  A(2)= CMPLX(-3.75d0,0.4d0)
  A(3)= CMPLX(4.25d0,-1.d0)

! Example #2 (real roots are -3, 1, 2)
! A(0)= CMPLX( 1.d0,0.d0)
! A(1)= CMPLX( 0.d0,0.d0)
! A(2)= CMPLX(-7.d0,0.d0)
! A(3)= CMPLX( 6.d0,0.d0)

  ITMAX=20                  !Maximum number of iterations
  EPS=1.d-8; EPS2=1.d-6     !Minimum, maximum relative error
  RAC(0)=CMPLX(0.d0,0.d0)   !Approximate value of 1st root

  call CLAGUE(N,A,ITMAX,EPS,EPS2,IMP,RAC,WORK)

  print *,' '
  print *,' Error code =', IMP
  print *,' '
  print *,' Complex roots are:'
  print *,' '
  do I=1, N
    write(*,*)  RAC(I)
  end do
  print *,' '

  stop

END
 
      SUBROUTINE CLAGUE(N,A,ITMAX,EPS,EPS2,IMP,X,W)
!=================================================================
!     ROOTS OF A COMPLEX COEFFICIENTS POLYNOMIAL
!     BY LAGUERRE FORMULA IN COMPLEX DOMAIN
!=================================================================
!     CALLING MODE:
!       CALL CLAGUE(N,A,ITMAX,EPS,EPS2,IMP,X,W)
!     INPUTS:
!     N : ORDER OF POLYNOMIAL
!     A : TABLE OF SIZE (0:N) OF COMPLEX COEFFICIENTS STORED IN
!         DECREASING ORDER OF X POWERS
!     ITMAX: MAXIMUN NUMBER OF ITERATIONS FOR EACH ROOT
!     EPS:   MINIMAL RELATIVE ERROR
!     EPS2:  MAXIMAL RELATIVE ERROR
!     X(0):  APPROXIMATE VALUE OF FIRST ROOT (=0. GENERALLY)
!     OUTPUTS:
!     IMP:  FLAG = 0  CONVERGENCE WITH AT LEAST EPS2 PRECISION
!                = 1  NO CONVERGENCE
!     X:  TABLE OF SIZE (0:N) OF FOUND ROOTS
!         X(I),I=1,N
!     WORKING ZONE:
!     W:  TABLE OF SIZE (0:3*N)
! ----------------------------------------------------------------
!     REFERENCE:
!     E.DURAND. SOLUTIONS NUMERIQUES DES EQUATIONS ALGEBRIQUES
!               TOME I, MASSON & CIE PAGES 269-270
!=================================================================
      COMPLEX A(0:N),X(0:N),W(0:3*N)
      REAL*8 EPS,EPS2
      N1=N+1
      N2=N1+N

      CALL LGUERR(N,A,ITMAX,EPS,EPS2,X,IMP,W(0),W(N1),W(N2))

      RETURN
      END

      SUBROUTINE LGUERR(K,A,ITMAX,EPS,EPS2,X,IMP,B,C,D)
      COMPLEX A(0:K),B(0:K),C(0:K),D(0:K),X(0:K),  &
      XK,XR,F,FP,FS,H,DEN,SQ,D2
      REAL*8 EPS,EPS2,TEST

      IMP=0
      B(0)=A(0)
      C(0)=B(0)
      D(0)=C(0)
      IF(K.EQ.1) THEN
      X(1)=-A(1)/A(0)
      RETURN
      ENDIF
      N=K

      XK=X(0)
    1 IK=0
      EPS1=EPS
      IT=0

    2 DO I=1,N
      B(I)=A(I)+XK*B(I-1)
      IF(I.LE.N-1) C(I)=B(I)+XK*C(I-1)
      IF(I.LE.N-2) D(I)=C(I)+XK*D(I-1)
      ENDDO

      F=B(N)
      FP=C(N-1)
      FS=D(N-2)
      H=(FLOAT(N-1)*FP)**2-FLOAT(N*(N-1))*F*FS
      IF(CABS(H).EQ.0.D0) THEN
      DO I=N,1,-1
      X(I)=-A(1)/(FLOAT(N)*A(0))
      ENDDO
      RETURN
      ENDIF
      IF(FP.EQ.0.D0) THEN
      DEN=CSQRT(H)
      ELSE
      SQ=CSQRT(H)
	  print *,' H=', H
	  print *,' SQ=', SQ
      DEN=FP+SQ
      D2=FP-SQ
      IF(CABS(DEN).LT.CABS(D2))  DEN=D2
      ENDIF
      IF(DEN.EQ.0.D0) THEN
      XK=XK+(0.1D0,0.1D0)
      GO TO 2
      ENDIF
      IK=IK+1
      XR=XK-FLOAT(N)*F/DEN
      print *,' XR-XK=', XR-XK
      TEST=CABS(XR-XK)
      XK=XR
      IF(TEST.LT.EPS1*CABS(XR)) GO TO 3
      WRITE(*,*) IK,H,DEN,XK,TEST
	  Pause
      IF(IK.LE.ITMAX)  GO TO 2
      EPS1=EPS1*10.D0
      IK=0
!     WRITE(*,*) 'EPS1=',EPS1
      IF(EPS1.LT.EPS2) GO TO 2
      IMP=1
      RETURN

    3 X(N)=XR
      N=N-1
      IF(N.EQ.1) GO TO 4
      IF(N.LE.0) RETURN
      DO I=0,N
      A(I)=B(I)
      ENDDO
      GO TO 1

    4 X(N)=-B(1)/B(0)
      RETURN
      END

!end of file tclague.f90