module defpi implicit none real(16),parameter::pi=4.0Q0*atan(1.0Q0) end module defpi program f90_DKA !*********************************************************** ! DKA method (Durand-Kerner-Aberth) ! f(x)=a0*x^n-a1*x^(n-1)+....+an=0 ! http://www.ecs.shimane-u.ac.jp/~kyoshida/c6(2002).pdf !*********************************************************** use defpi implicit none integer,parameter::itmax=10000 real(16),parameter::eps=1.0Q-20 real(16),parameter::zero=1.0Q-1000 integer::i,j,n,iii real(16),allocatable::a(:),z1(:),z2(:),err(:) real(16),allocatable::ff0(:),z1mem(:),z2mem(:) integer,allocatable::imem(:),icoef(:) real(16)::w0,w1,w2,f1,f2,ffe real(16)::p1,p2,ee,ae,a1,a2,sume real(16)::rr,rc,theta,xc,yc character(len=50)::fnameR character(len=100)::str !---------------------------------------------------- ! Data input !---------------------------------------------------- call getarg(1,fnameR) open(11,file=fnameR,status='old') read(11,'(a)') str read(11,*) n allocate(a(0:n),z1(1:n),z2(1:n),err(1:n)) allocate(ff0(1:n),z1mem(1:n),z2mem(1:n)) allocate(imem(1:n),icoef(1:n)) do i=0,n read(11,*) a(i) end do close(11) !---------------------------------------------------- ! Treatment of coefficient !---------------------------------------------------- w0=a(0) do i=0,n a(i)=a(i)/w0 end do !---------------------------------------------------- ! Initial value of x !---------------------------------------------------- write(6,'(a)') trim(adjustl(str)) ! rc=RINI1(n,a) rc=RINI2(n,a) do i=1,n theta=2.0Q0*pi/dble(n)*(dble(i)-0.75Q0) z1(i)=-a(1)/dble(n)+rc*cos(theta) z2(i)=rc*sin(theta) end do !---------------------------------------------------- ! Output of initial value !---------------------------------------------------- ! Write to the standard output write(6,'(a)') 'Initial value' write(6,'(a)') ' rc -a(1)/n' write(6,'(2(e15.7))') rc,-a(1)/dble(n) write(6,'(a)') ' i z1 z2 rr' do i=1,n rr=sqrt((z1(i)+a(1)/dble(n))**2+z2(i)**2) write(6,'(i5,3(e15.7))') i,z1(i),z2(i),rr end do ! Write to the file for GMT open(12, file='_iteration.txt', status='replace') do i=1,n write(12,'(2(e15.7))') z1(i),z2(i) end do !---------------------------------------------------- ! Iteration by DK method !---------------------------------------------------- do i=1,n ff0(i)=1.0Q1000 end do do iii=1,itmax ee=0.0Q0 sume=0.0Q0 do i=1,n ! Numerator f1=1.0Q0 f2=0.0Q0 do j=1,n w1=f1*z1(i)-f2*z2(i) w2=f2*z1(i)+f1*z2(i) f1=w1+a(j) f2=w2 end do ! Denominator p1=1.0Q0 p2=0.0Q0 do j=1,n if(j==i)cycle w1=p1*(z1(i)-z1(j))-p2*(z2(i)-z2(j)) w2=p1*(z2(i)-z2(j))+p2*(z1(i)-z1(j)) p1=w1 p2=w2 end do a1=(f1*p1+f2*p2)/(p1*p1+p2*p2) a2=(f2*p1-f1*p2)/(p1*p1+p2*p2) ! Evaluation of (f/df) ae=sqrt(a1*a1+a2*a2) err(i)=ae if(ee')) do k=0,iii write(12,'(2(e15.7))') zz1(k,i),zz2(k,i) end do end do close(12) open(12, file='_data_circ.txt', status='replace') do i=0,360 zx=zxc+rr*cos(dble(i)/180.0Q0*pi) zy=zyc+rr*sin(dble(i)/180.0Q0*pi) write(12,'(2(e15.7))') zx,zy end do close(12) end subroutine GMTDAT !---------------------------------------------------- ! Calculation of r-value for initial x-value ! (simple method) !---------------------------------------------------- real(16) function RINI1(n,a) integer,intent(in)::n real(16),intent(in)::a(0:n) integer::i,ic real(16)::rm ic=0 do i=1,n if(eps