module fvn_misc use fvn_common implicit none ! Muller interface fvn_muller module procedure fvn_z_muller end interface fvn_muller contains ! ! Muller ! ! ! ! William Daniau 2007 ! ! This routine is a fortran 90 port of Hans D. Mittelmann's routine muller.f ! http://plato.asu.edu/ftp/other_software/muller.f ! ! it can be used as a replacement for imsl routine dzanly with minor changes ! !----------------------------------------------------------------------- ! ! purpose - zeros of an analytic complex function ! using the muller method with deflation ! ! usage - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax, ! infer,ier) ! ! arguments f - a complex function subprogram, f(z), written ! by the user specifying the equation whose ! roots are to be found. f must appear in ! an external statement in the calling pro- ! gram. ! eps - 1st stopping criterion. let fp(z)=f(z)/p ! where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1)) ! and z(1),...,z(k-1) are previously found ! roots. if ((cdabs(f(z)).le.eps) .and. ! (cdabs(fp(z)).le.eps)), then z is accepted ! as a root. (input) ! eps1 - 2nd stopping criterion. a root is accepted ! if two successive approximations to a given ! root agree within eps1. (input) ! note. if either or both of the stopping ! criteria are fulfilled, the root is ! accepted. ! kn - the number of known roots which must be stored ! in x(1),...,x(kn), prior to entry to muller ! nguess - the number of initial guesses provided. these ! guesses must be stored in x(kn+1),..., ! x(kn+nguess). nguess must be set equal ! to zero if no guesses are provided. (input) ! n - the number of new roots to be found by ! muller (input) ! x - a complex vector of length kn+n. x(1),..., ! x(kn) on input must contain any known ! roots. x(kn+1),..., x(kn+n) on input may, ! on user option, contain initial guesses for ! the n new roots which are to be computed. ! if the user does not provide an initial ! guess, zero is used. ! on output, x(kn+1),...,x(kn+n) contain the ! approximate roots found by muller. ! itmax - the maximum allowable number of iterations ! per root (input) ! infer - an integer vector of length kn+n. on ! output infer(j) contains the number of ! iterations used in finding the j-th root ! when convergence was achieved. if ! convergence was not obtained in itmax ! iterations, infer(j) will be greater than ! itmax (output). ! ier - error parameter (output) ! warning error ! ier = 33 indicates failure to converge with- ! in itmax iterations for at least one of ! the (n) new roots. ! ! ! remarks muller always returns the last approximation for root j ! in x(j). if the convergence criterion is satisfied, ! then infer(j) is less than or equal to itmax. if the ! convergence criterion is not satisified, then infer(j) ! is set to either itmax+1 or itmax+k, with k greater ! than 1. infer(j) = itmax+1 indicates that muller did ! not obtain convergence in the allowed number of iter- ! ations. in this case, the user may wish to set itmax ! to a larger value. infer(j) = itmax+k means that con- ! vergence was obtained (on iteration k) for the defla- ! ted function ! fp(z) = f(z)/((z-z(1)...(z-z(j-1))) ! ! but failed for f(z). in this case, better initial ! guesses might help or, it might be necessary to relax ! the convergence criterion. ! !----------------------------------------------------------------------- ! subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier) implicit none double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq,eps1w double complex :: d,dd,den,fprt,frt,h,rt,t1,t2,t3, & tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, & zero,p1,one,four,p5 double complex, external :: f integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, & knpng,jk,ick,nn,lm1,errcode double complex :: x(kn+n) integer :: infer(kn+n) data zero/(0.0,0.0)/,p1/(0.1,0.0)/, & one/(1.0,0.0)/,four/(4.0,0.0)/, & p5/(0.5,0.0)/, & rzero/0.0/,rten/10.0/,rhun/100.0/, & ax/0.1/,ickmax/3/,rp01/0.01/ ier = 0 if (n .lt. 1) then ! What the hell are doing here then ... return end if !eps1 = rten **(-nsig) eps1w = min(eps1,rp01) knp1 = kn+1 knpn = kn+n knpng = kn+nguess do i=1,knpn infer(i) = 0 if (i .gt. knpng) x(i) = zero end do l= knp1 ic=0 rloop: do while (l<=knpn) ! Main loop over new roots jk = 0 ick = 0 xl = x(l) icloop: do ic = 0 h = ax h = p1*h if (cdabs(xl) .gt. ax) h = p1*xl ! first three points are ! xl+h, xl-h, xl rt = xl+h call deflated_work(errcode) if (errcode == 1) then exit icloop end if z0 = fprt y0 = frt x0 = rt rt = xl-h call deflated_work(errcode) if (errcode == 1) then exit icloop end if z1 = fprt y1 = frt h = xl-rt d = h/(rt-x0) rt = xl call deflated_work(errcode) if (errcode == 1) then exit icloop end if z2 = fprt y2 = frt ! begin main algorithm iloop: do dd = one + d t1 = z0*d*d t2 = z1*dd*dd xx = z2*dd t3 = z2*d bi = t1-t2+xx+t3 den = bi*bi-four*(xx*t1-t3*(t2-xx)) ! use denominator of maximum amplitude t1 = cdsqrt(den) qz = rhun*max(cdabs(bi),cdabs(t1)) t2 = bi + t1 tpq = cdabs(t2)+qz if (tpq .eq. qz) t2 = zero t3 = bi - t1 tpq = cdabs(t3) + qz if (tpq .eq. qz) t3 = zero den = t2 qz = cdabs(t3)-cdabs(t2) if (qz .gt. rzero) den = t3 ! test for zero denominator if (cdabs(den) .eq. rzero) then call trans_rt() call deflated_work(errcode) if (errcode == 1) then exit icloop end if z2 = fprt y2 = frt cycle iloop end if d = -xx/den d = d+d h = d*h rt = rt + h ! check convergence of the first kind if (cdabs(h) .le. eps1w*max(cdabs(rt),ax)) then if (ic .ne. 0) then exit icloop end if ic = 1 z0 = y1 z1 = y2 z2 = f(rt) xl = rt ick = ick+1 if (ick .le. ickmax) then cycle iloop end if ! warning error, itmax = maximum jk = itmax + jk ier = 33 end if if (ic .ne. 0) then cycle icloop end if call deflated_work(errcode) if (errcode == 1) then exit icloop end if do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero) ! take remedial action to induce ! convergence d = d*p5 h = h*p5 rt = rt-h call deflated_work(errcode) if (errcode == 1) then exit icloop end if end do z0 = z1 z1 = z2 z2 = fprt y0 = y1 y1 = y2 y2 = frt end do iloop end do icloop x(l) = rt infer(l) = jk l = l+1 end do rloop contains subroutine trans_rt() tem = rten*eps1w if (cdabs(rt) .gt. ax) tem = tem*rt rt = rt+tem d = (h+tem)*d/h h = h+tem end subroutine trans_rt subroutine deflated_work(errcode) ! errcode=0 => no errors ! errcode=1 => jk>itmax or convergence of second kind achieved integer :: errcode,flag flag=1 loop1: do while(flag==1) errcode=0 jk = jk+1 if (jk .gt. itmax) then ier=33 errcode=1 return end if frt = f(rt) fprt = frt if (l /= 1) then lm1 = l-1 do i=1,lm1 tem = rt - x(i) if (cdabs(tem) .eq. rzero) then !if (ic .ne. 0) go to 15 !! ?? possible? call trans_rt() cycle loop1 end if fprt = fprt/tem end do end if flag=0 end do loop1 if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then errcode=1 return end if end subroutine deflated_work end subroutine end module fvn_misc