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.0d0,0.0d0)/,p1/(0.1d0,0.0d0)/, &
                          one/(1.0d0,0.0d0)/,four/(4.0d0,0.0d0)/, &
                          p5/(0.5d0,0.0d0)/, &
                          rzero/0.0d0/,rten/10.0d0/,rhun/100.0d0/, &
                          ax/0.1d0/,ickmax/3/,rp01/0.01d0/

            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 = sqrt(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