source/unres/quindibisectok.F90

   1         subroutine quindibisect2(c,b,d,n,m1,m2,eps1,relfeh,eps2,z,x)
   2         real (kind=8) :: eps1,eps2,relfeh,xmin,xmax,h,x1,xu,xo
   3         real (kind=8),dimension(n):: c,b,d,x
   4         real (kind=8) ,dimension(:):: wu(m1:m2)
   5         integer :: i,n,m1,m2,z,a,k
   6         integer maxit
   7         parameter (maxit=1000)
   8
   9         d(1)=0
  10         d(2)=0
  11         b(1)=0
  12         xmin = c(n)-dabs(b(n))-dabs(d(n))
  13         xmax = c(n) + dabs(b(n)) + dabs(d(n))
  14         h=dabs(b(n-1))+dabs(d(n-1))+dabs(b(n))
  15         if ((c(n-1)+h).gt.xmax) xmax=c(n-1)+h
  16         if ((c(n-1)-h).lt.xmin) xmin=c(n-1)-h
  17         do i=n-2,1,-1
  18             h= dabs(b(i))+dabs(d(i))+dabs(b(i+1))+dabs(d(i+2))
  19             if ((c(i)+h).gt.xmax) xmax=c(i)+h
  20             if ((c(i)-h).lt.xmin) xmin=c(i)-h
  21         enddo
  22         if ((xmin+xmax).gt.0) then
  23            eps2=relfeh*xmax
  24         else
  25            eps2=-relfeh*xmin
  26         endif
  27         if (eps1.lt.0) eps1=eps2
  28         eps2= 0.5*eps1+7*eps2
  29         xo=xmax
  30         do i=m1,m2
  31           x(i)=xmax
  32           wu(i)=xmin
  33         enddo
  34         z=0
  35         do k=m2,m1,-1
  36           xu=xmin
  37 !comment loop for the kth value
  38           do i=k,m1,-1
  39             if (xu.lt.wu(i)) then
  40               xu=wu(i)
  41               go to 10
  42             endif
  43           enddo
  44        10 if (xo.gt.x(k)) xo=x(k)
  45           x1=(xu+xo)/2
  46           do while  (xo-xu.gt.2*relfeh*(dabs(xu)+dabs(xo)+eps1))
  47              if (z.gt.maxit) exit
  48              z=z+1
  49              call quindet2(c,b,d,n,x1,a)
  50              a=n-a
  51              if (a.lt.k) then
  52                 if (a.lt.m1) then
  53                    wu(m1)=x1
  54                    xu=x1
  55                 else
  56                   wu(a+1)=x1
  57                   xu=x1
  58                   if (x(a).gt.x1) x(a)=x1
  59                 endif
  60              else
  61                 xo=x1
  62              endif
  63              x1=(xu+xo)/2
  64           enddo
  65           x(k)=(xo+xu)/2
  66         enddo ! k
  67        endsubroutine
  68