program rootf * * This is a polynomial root finder program. * The polynomials must have real coefficients. * The method applied is a modification of BAIRSTOW's * i.e. division by a quadratic polynomial x**2 +p1*x +p2 * determining p1 & p2 to be such that the remainder is vanishing. * Then a deflation is applied by synthetic division to reduce * the polynomial degree by two, and the procedure is repeated. * A check is always made if the degree is reduced to either * one or two to use special finishing procedure * IMPLICIT DOUBLE PRECISION (A-H,O-Z) common/toroot/a(0:100),b(0:100),kco COMPLEX ROOT(100),disc parameter (N=2,M=2) dimension icode(4),xp(2),xll(2),xrl(2),ixat(2) character line*78 data icode/1,3*0/ do 100 i=1,78 100 line(I:I)='=' * Read in, tolerance, polynomial degree, polynomial coefficients write(*,*) 'Enter a termination tolerance: ' read(*,*) eps write(*,*) 'Enter the degree of the polynomial:' read(*,*)kco write(*,*) 'Enter the polynomial coefficients:( a(0),a(1),...)' read(*,*)(a(i),i=0,kco) * Type the above input for reference write(*,'(A)')line write(*,*) ' RUNNING WITH THE FOLLOWING INPUT ' write(*,*) 'TOLERANCE: ',eps write(*,*) 'DEGREE: ',kco write(*,*)' COEFFICIENTS a(i), i=0,1,... ' write(*,*)(a(i),i=0,kco) write(*,'(A)')line nor = 0 30 continue * Check if the degree is one, and if so solve a linear equation. if (kco.eq.1) then write(*,*)' The roots of the polynomial are:' nor = nor + 1 root(nor)= cmplx(-a(0)/a(1),0) do 1 i=1,nor write(*,19)i,root(i) 1 continue 19 format(2x,'Root # ',i4,t17,'(',g14.7,', ',g14.7,')') stop endif * Check if the degree is two, and if so solve a quadratic equation. if (kco .eq. 2) then write(*,*)' The roots of the polynomial are:' nor = nor + 1 disc = cmplx(a(1)**2 - 4*a(2)*a(0),0.d0) root(nor) = (-a(1)+csqrt(disc))/2./a(2) nor = nor + 1 root(nor) = (-a(1)-csqrt(disc))/2./a(2) do 2 i=1,nor write(*,19)i,root(i) 2 continue stop endif * Initialize randomly in (-1,1) p1 = 2*ranm()-1 p2 = 2*ranm()-1 xp(1) = p1 xp(2) = p2 * Minimize the remainder down to zero !!! call OPTIMA(N,M,XP,XV,XLL,XRL,IXAT,ICODE, > 'in.dat','out.dat',grms,NF,NG,NHE,NJA) if (sqrt(xv) .le. eps) then p1 = xp(1) p2 = xp(2) disc = cmplx( p1**2-4*p2, 0.d0) nor = nor + 1 root(nor)=(-p1 + csqrt(disc))/2 nor = nor + 1 root(nor)=(-p1 - csqrt(disc))/2 call syndiv(a,kco,p1,p2,b,R,Q) kco = kco - 2 do 3 i = 0,kco a(i) = b(i) 3 continue endif go to 30 end * subroutine syndiv(a,n,p1,p2,b,R,Q) IMPLICIT DOUBLE PRECISION (A-H,O-Z) * Given a Polynomial: a(n)*x**n + ...+a(1)*x +a(0) * and a quadratic: x**2 + p1*x + p2 * calculate the coefficients of the quotient polynomial: b(k),k=0,1,...,n-2 * and the remainder terms R and Q. dimension a(0:*),b(0:*) b(n-2)=a(n) b(n-3)=a(n-1)-p1*b(n-2) do 1 k = n-2,2,-1 b(k-2) = a(k)-p1*b(k-1)-p2*b(k) 1 continue R = a(1) - p1*b(0)-p2*b(1) Q = a(0) - p2*b(0) end * subroutine synder(a,n,p1,p2,b,R,Q,b1,b2,R1,R2,Q1,Q2) IMPLICIT DOUBLE PRECISION (A-H,O-Z) * Given a Polynomial: a(n)*x**n + ...+a(1)*x +a(0) * and a quadratic: x**2 + p1*x + p2 * Calculate the coefficients of the quotient polynomial: b(k),k=0,1,...,n-2 * and the remainder terms R and Q. * Calculate the derivatives: b1(k),R1,Q1 of b(k),R and Q with respect to p1, * and the derivatives: b2(k),R2,Q2 of b(k),R and Q with respect to p2. dimension a(0:*),b(0:*),b1(0:*),b2(0:*) * b(n-2)=a(n) b(n-3)=a(n-1) - p1*b(n-2) do 1 k = n-2,2,-1 b(k-2) = a(k) - p1*b(k-1) - p2*b(k) 1 continue R = a(1) - p1*b(0)-p2*b(1) Q = a(0) - p2*b(0) * b1(n-2)= 0 b1(n-3)= -a(n) do 2 k = n-2,2,-1 b1(k-2) = -b(k-1) - p1*b1(k-1) - p2*b1(k) 2 continue R1 = -b(0) - p1*b1(0) - p2*b1(1) Q1 = -p2*b1(0) * b2(n-2)= 0 b2(n-3)= 0 do 3 k = n-2,2,-1 b2(k-2) = -p1*b2(k-1) - b(k) - p2*b2(k) 3 continue R2 = -p1*b2(0) - b(1) - p2*b2(1) Q2 = -b(0) - p2*b2(0) end C --------------------------------------------------------------------- SUBROUTINE SUBSUM ( M, N, X, F ) C --------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(N), F(M) common/toroot/a(0:100),b(0:100),kco p1 = x(1) p2 = x(2) call syndiv(a,kco,p1,p2,b,R,Q) f(1) = R f(2) = Q END C --------------------------------------------------------------------- SUBROUTINE GRANAL ( N, X, GRAD ) C --------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(N), GRAD(N) common/toroot/a(0:100),b(0:100),kco dimension b1(0:100),b2(0:100) p1 = x(1) p2 = x(2) call synder(a,kco,p1,p2,b,R,Q,b1,b2,R1,R2,Q1,Q2) grad(1) = 2*(R*R1 + Q*Q1) grad(2) = 2*(R*R2 + Q*Q2) end C --------------------------------------------------------------------- SUBROUTINE JANAL ( M, N, X, FJ, LD ) C --------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(N), FJ(LD,N) common/toroot/a(0:100),b(0:100),kco dimension b1(0:100),b2(0:100) p1 = x(1) p2 = x(2) call synder(a,kco,p1,p2,b,R,Q,b1,b2,R1,R2,Q1,Q2) FJ(1,1) = R1 FJ(1,2) = R2 FJ(2,1) = Q1 FJ(2,2) = Q2 END