Blame view
fvn_misc/fvn_misc.f90
12 KB
b93026039
git-svn-id: https...
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 |
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
|