!**************************************************** !* 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