lxsd / fvn

Download as
Browse Code ยป

Commit 9c285563c30716b63de00e330bdfe930a39a21a3

Authored by kwagner
1 parent 00055ac083
Exists in master and in 2 other branches geevx_2020, multi

Define constants in data statements as double precision 

git-svn-id: https://lxsd.femto-st.fr/svn/fvn@71 b657c933-2333-4658-acf2-d3c7c2708721

Showing 1 changed file with 7 additions and 7 deletions Inline Diff

fvn_misc/fvn_misc.f90
Diff comments   View file @ 9c28556
module fvn_misc 1 1 module fvn_misc
use fvn_common 2 2 use fvn_common
implicit none 3 3 implicit none
4 4
! Muller 5 5 ! Muller
interface fvn_muller 6 6 interface fvn_muller
module procedure fvn_z_muller 7 7 module procedure fvn_z_muller
end interface fvn_muller 8 8 end interface fvn_muller
9 9
contains 10 10 contains
! 11 11 !
! Muller 12 12 ! Muller
! 13 13 !
! 14 14 !
! 15 15 !
! William Daniau 2007 16 16 ! William Daniau 2007
! 17 17 !
! This routine is a fortran 90 port of Hans D. Mittelmann's routine muller.f 18 18 ! This routine is a fortran 90 port of Hans D. Mittelmann's routine muller.f
! http://plato.asu.edu/ftp/other_software/muller.f 19 19 ! http://plato.asu.edu/ftp/other_software/muller.f
! 20 20 !
! it can be used as a replacement for imsl routine dzanly with minor changes 21 21 ! it can be used as a replacement for imsl routine dzanly with minor changes
! 22 22 !
!----------------------------------------------------------------------- 23 23 !-----------------------------------------------------------------------
! 24 24 !
! purpose - zeros of an analytic complex function 25 25 ! purpose - zeros of an analytic complex function
! using the muller method with deflation 26 26 ! using the muller method with deflation
! 27 27 !
! usage - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax, 28 28 ! usage - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax,
! infer,ier) 29 29 ! infer,ier)
! 30 30 !
! arguments f - a complex function subprogram, f(z), written 31 31 ! arguments f - a complex function subprogram, f(z), written
! by the user specifying the equation whose 32 32 ! by the user specifying the equation whose
! roots are to be found. f must appear in 33 33 ! roots are to be found. f must appear in
! an external statement in the calling pro- 34 34 ! an external statement in the calling pro-
! gram. 35 35 ! gram.
! eps - 1st stopping criterion. let fp(z)=f(z)/p 36 36 ! eps - 1st stopping criterion. let fp(z)=f(z)/p
! where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1)) 37 37 ! where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1))
! and z(1),...,z(k-1) are previously found 38 38 ! and z(1),...,z(k-1) are previously found
! roots. if ((cdabs(f(z)).le.eps) .and. 39 39 ! roots. if ((cdabs(f(z)).le.eps) .and.
! (cdabs(fp(z)).le.eps)), then z is accepted 40 40 ! (cdabs(fp(z)).le.eps)), then z is accepted
! as a root. (input) 41 41 ! as a root. (input)
! eps1 - 2nd stopping criterion. a root is accepted 42 42 ! eps1 - 2nd stopping criterion. a root is accepted
! if two successive approximations to a given 43 43 ! if two successive approximations to a given
! root agree within eps1. (input) 44 44 ! root agree within eps1. (input)
! note. if either or both of the stopping 45 45 ! note. if either or both of the stopping
! criteria are fulfilled, the root is 46 46 ! criteria are fulfilled, the root is
! accepted. 47 47 ! accepted.
! kn - the number of known roots which must be stored 48 48 ! kn - the number of known roots which must be stored
! in x(1),...,x(kn), prior to entry to muller 49 49 ! in x(1),...,x(kn), prior to entry to muller
! nguess - the number of initial guesses provided. these 50 50 ! nguess - the number of initial guesses provided. these
! guesses must be stored in x(kn+1),..., 51 51 ! guesses must be stored in x(kn+1),...,
! x(kn+nguess). nguess must be set equal 52 52 ! x(kn+nguess). nguess must be set equal
! to zero if no guesses are provided. (input) 53 53 ! to zero if no guesses are provided. (input)
! n - the number of new roots to be found by 54 54 ! n - the number of new roots to be found by
! muller (input) 55 55 ! muller (input)
! x - a complex vector of length kn+n. x(1),..., 56 56 ! x - a complex vector of length kn+n. x(1),...,
! x(kn) on input must contain any known 57 57 ! x(kn) on input must contain any known
! roots. x(kn+1),..., x(kn+n) on input may, 58 58 ! roots. x(kn+1),..., x(kn+n) on input may,
! on user option, contain initial guesses for 59 59 ! on user option, contain initial guesses for
! the n new roots which are to be computed. 60 60 ! the n new roots which are to be computed.
! if the user does not provide an initial 61 61 ! if the user does not provide an initial
! guess, zero is used. 62 62 ! guess, zero is used.
! on output, x(kn+1),...,x(kn+n) contain the 63 63 ! on output, x(kn+1),...,x(kn+n) contain the
! approximate roots found by muller. 64 64 ! approximate roots found by muller.
! itmax - the maximum allowable number of iterations 65 65 ! itmax - the maximum allowable number of iterations
! per root (input) 66 66 ! per root (input)
! infer - an integer vector of length kn+n. on 67 67 ! infer - an integer vector of length kn+n. on
! output infer(j) contains the number of 68 68 ! output infer(j) contains the number of
! iterations used in finding the j-th root 69 69 ! iterations used in finding the j-th root
! when convergence was achieved. if 70 70 ! when convergence was achieved. if
! convergence was not obtained in itmax 71 71 ! convergence was not obtained in itmax
! iterations, infer(j) will be greater than 72 72 ! iterations, infer(j) will be greater than
! itmax (output). 73 73 ! itmax (output).
! ier - error parameter (output) 74 74 ! ier - error parameter (output)
! warning error 75 75 ! warning error
! ier = 33 indicates failure to converge with- 76 76 ! ier = 33 indicates failure to converge with-
! in itmax iterations for at least one of 77 77 ! in itmax iterations for at least one of
! the (n) new roots. 78 78 ! the (n) new roots.
! 79 79 !
! 80 80 !
! remarks muller always returns the last approximation for root j 81 81 ! remarks muller always returns the last approximation for root j
! in x(j). if the convergence criterion is satisfied, 82 82 ! in x(j). if the convergence criterion is satisfied,
! then infer(j) is less than or equal to itmax. if the 83 83 ! then infer(j) is less than or equal to itmax. if the
! convergence criterion is not satisified, then infer(j) 84 84 ! convergence criterion is not satisified, then infer(j)
! is set to either itmax+1 or itmax+k, with k greater 85 85 ! is set to either itmax+1 or itmax+k, with k greater
! than 1. infer(j) = itmax+1 indicates that muller did 86 86 ! than 1. infer(j) = itmax+1 indicates that muller did
! not obtain convergence in the allowed number of iter- 87 87 ! not obtain convergence in the allowed number of iter-
! ations. in this case, the user may wish to set itmax 88 88 ! ations. in this case, the user may wish to set itmax
! to a larger value. infer(j) = itmax+k means that con- 89 89 ! to a larger value. infer(j) = itmax+k means that con-
! vergence was obtained (on iteration k) for the defla- 90 90 ! vergence was obtained (on iteration k) for the defla-
! ted function 91 91 ! ted function
! fp(z) = f(z)/((z-z(1)...(z-z(j-1))) 92 92 ! fp(z) = f(z)/((z-z(1)...(z-z(j-1)))
! 93 93 !
! but failed for f(z). in this case, better initial 94 94 ! but failed for f(z). in this case, better initial
! guesses might help or, it might be necessary to relax 95 95 ! guesses might help or, it might be necessary to relax
! the convergence criterion. 96 96 ! the convergence criterion.
! 97 97 !
!----------------------------------------------------------------------- 98 98 !-----------------------------------------------------------------------
! 99 99 !
subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier) 100 100 subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)
implicit none 101 101 implicit none
double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq,eps1w 102 102 double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq,eps1w
double complex :: d,dd,den,fprt,frt,h,rt,t1,t2,t3, & 103 103 double complex :: d,dd,den,fprt,frt,h,rt,t1,t2,t3, &
tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, & 104 104 tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, &
zero,p1,one,four,p5 105 105 zero,p1,one,four,p5
106 106
double complex, external :: f 107 107 double complex, external :: f
integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, & 108 108 integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, &
knpng,jk,ick,nn,lm1,errcode 109 109 knpng,jk,ick,nn,lm1,errcode
double complex :: x(kn+n) 110 110 double complex :: x(kn+n)
integer :: infer(kn+n) 111 111 integer :: infer(kn+n)
112 112
113 113
data zero/(0.0,0.0)/,p1/(0.1,0.0)/, & 114 114 data zero/(0.0d0,0.0d0)/,p1/(0.1d0,0.0d0)/, &
one/(1.0,0.0)/,four/(4.0,0.0)/, & 115 115 one/(1.0d0,0.0d0)/,four/(4.0d0,0.0d0)/, &
p5/(0.5,0.0)/, & 116 116 p5/(0.5d0,0.0d0)/, &
rzero/0.0/,rten/10.0/,rhun/100.0/, & 117 117 rzero/0.0d0/,rten/10.0d0/,rhun/100.0d0/, &
ax/0.1/,ickmax/3/,rp01/0.01/ 118 118 ax/0.1d0/,ickmax/3/,rp01/0.01d0/
119 119
ier = 0 120 120 ier = 0
if (n .lt. 1) then ! What the hell are doing here then ... 121 121 if (n .lt. 1) then ! What the hell are doing here then ...
return 122 122 return
end if 123 123 end if
!eps1 = rten **(-nsig) 124 124 !eps1 = rten **(-nsig)
eps1w = min(eps1,rp01) 125 125 eps1w = min(eps1,rp01)
126 126
knp1 = kn+1 127 127 knp1 = kn+1
knpn = kn+n 128 128 knpn = kn+n
knpng = kn+nguess 129 129 knpng = kn+nguess
do i=1,knpn 130 130 do i=1,knpn
infer(i) = 0 131 131 infer(i) = 0
if (i .gt. knpng) x(i) = zero 132 132 if (i .gt. knpng) x(i) = zero
end do 133 133 end do
l= knp1 134 134 l= knp1
135 135
ic=0 136 136 ic=0
rloop: do while (l<=knpn) ! Main loop over new roots 137 137 rloop: do while (l<=knpn) ! Main loop over new roots
jk = 0 138 138 jk = 0
ick = 0 139 139 ick = 0
xl = x(l) 140 140 xl = x(l)
icloop: do 141 141 icloop: do
ic = 0 142 142 ic = 0
h = ax 143 143 h = ax
h = p1*h 144 144 h = p1*h
if (cdabs(xl) .gt. ax) h = p1*xl 145 145 if (cdabs(xl) .gt. ax) h = p1*xl
! first three points are 146 146 ! first three points are
! xl+h, xl-h, xl 147 147 ! xl+h, xl-h, xl
rt = xl+h 148 148 rt = xl+h
call deflated_work(errcode) 149 149 call deflated_work(errcode)
if (errcode == 1) then 150 150 if (errcode == 1) then
exit icloop 151 151 exit icloop
end if 152 152 end if
153 153
z0 = fprt 154 154 z0 = fprt
y0 = frt 155 155 y0 = frt
x0 = rt 156 156 x0 = rt
rt = xl-h 157 157 rt = xl-h
call deflated_work(errcode) 158 158 call deflated_work(errcode)
if (errcode == 1) then 159 159 if (errcode == 1) then
exit icloop 160 160 exit icloop
end if 161 161 end if
162 162
z1 = fprt 163 163 z1 = fprt
y1 = frt 164 164 y1 = frt
h = xl-rt 165 165 h = xl-rt
d = h/(rt-x0) 166 166 d = h/(rt-x0)
rt = xl 167 167 rt = xl
168 168
call deflated_work(errcode) 169 169 call deflated_work(errcode)
if (errcode == 1) then 170 170 if (errcode == 1) then
exit icloop 171 171 exit icloop
end if 172 172 end if
173 173
174 174
z2 = fprt 175 175 z2 = fprt
y2 = frt 176 176 y2 = frt
! begin main algorithm 177 177 ! begin main algorithm
iloop: do 178 178 iloop: do
dd = one + d 179 179 dd = one + d
t1 = z0*d*d 180 180 t1 = z0*d*d
t2 = z1*dd*dd 181 181 t2 = z1*dd*dd
xx = z2*dd 182 182 xx = z2*dd
t3 = z2*d 183 183 t3 = z2*d
bi = t1-t2+xx+t3 184 184 bi = t1-t2+xx+t3
den = bi*bi-four*(xx*t1-t3*(t2-xx)) 185 185 den = bi*bi-four*(xx*t1-t3*(t2-xx))
! use denominator of maximum amplitude 186 186 ! use denominator of maximum amplitude
t1 = cdsqrt(den) 187 187 t1 = sqrt(den)
qz = rhun*max(cdabs(bi),cdabs(t1)) 188 188 qz = rhun*max(cdabs(bi),cdabs(t1))
t2 = bi + t1 189 189 t2 = bi + t1
tpq = cdabs(t2)+qz 190 190 tpq = cdabs(t2)+qz
if (tpq .eq. qz) t2 = zero 191 191 if (tpq .eq. qz) t2 = zero
t3 = bi - t1 192 192 t3 = bi - t1
tpq = cdabs(t3) + qz 193 193 tpq = cdabs(t3) + qz
if (tpq .eq. qz) t3 = zero 194 194 if (tpq .eq. qz) t3 = zero
den = t2 195 195 den = t2
qz = cdabs(t3)-cdabs(t2) 196 196 qz = cdabs(t3)-cdabs(t2)
if (qz .gt. rzero) den = t3 197 197 if (qz .gt. rzero) den = t3
! test for zero denominator 198 198 ! test for zero denominator
if (cdabs(den) .eq. rzero) then 199 199 if (cdabs(den) .eq. rzero) then
call trans_rt() 200 200 call trans_rt()
call deflated_work(errcode) 201 201 call deflated_work(errcode)
if (errcode == 1) then 202 202 if (errcode == 1) then
exit icloop 203 203 exit icloop
end if 204 204 end if
z2 = fprt 205 205 z2 = fprt
y2 = frt 206 206 y2 = frt
cycle iloop 207 207 cycle iloop
end if 208 208 end if
209 209
210 210
d = -xx/den 211 211 d = -xx/den
d = d+d 212 212 d = d+d
h = d*h 213 213 h = d*h
rt = rt + h 214 214 rt = rt + h
! check convergence of the first kind 215 215 ! check convergence of the first kind
if (cdabs(h) .le. eps1w*max(cdabs(rt),ax)) then 216 216 if (cdabs(h) .le. eps1w*max(cdabs(rt),ax)) then
if (ic .ne. 0) then 217 217 if (ic .ne. 0) then
exit icloop 218 218 exit icloop
end if 219 219 end if
ic = 1 220 220 ic = 1
z0 = y1 221 221 z0 = y1
z1 = y2 222 222 z1 = y2
z2 = f(rt) 223 223 z2 = f(rt)
xl = rt 224 224 xl = rt
ick = ick+1 225 225 ick = ick+1
if (ick .le. ickmax) then 226 226 if (ick .le. ickmax) then
cycle iloop 227 227 cycle iloop
end if 228 228 end if
! warning error, itmax = maximum 229 229 ! warning error, itmax = maximum
jk = itmax + jk 230 230 jk = itmax + jk
ier = 33 231 231 ier = 33
end if 232 232 end if
if (ic .ne. 0) then 233 233 if (ic .ne. 0) then
cycle icloop 234 234 cycle icloop
end if 235 235 end if
call deflated_work(errcode) 236 236 call deflated_work(errcode)
if (errcode == 1) then 237 237 if (errcode == 1) then
exit icloop 238 238 exit icloop
end if 239 239 end if
240 240
do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero) 241 241 do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero)
! take remedial action to induce 242 242 ! take remedial action to induce
! convergence 243 243 ! convergence
d = d*p5 244 244 d = d*p5
h = h*p5 245 245 h = h*p5
rt = rt-h 246 246 rt = rt-h
call deflated_work(errcode) 247 247 call deflated_work(errcode)
if (errcode == 1) then 248 248 if (errcode == 1) then
exit icloop 249 249 exit icloop
end if 250 250 end if
end do 251 251 end do
z0 = z1 252 252 z0 = z1
z1 = z2 253 253 z1 = z2
z2 = fprt 254 254 z2 = fprt
y0 = y1 255 255 y0 = y1
y1 = y2 256 256 y1 = y2
y2 = frt 257 257 y2 = frt
end do iloop 258 258 end do iloop
end do icloop 259 259 end do icloop
x(l) = rt 260 260 x(l) = rt
infer(l) = jk 261 261 infer(l) = jk
l = l+1 262 262 l = l+1
end do rloop 263 263 end do rloop
264 264
contains 265 265 contains
subroutine trans_rt() 266 266 subroutine trans_rt()
tem = rten*eps1w 267 267 tem = rten*eps1w
if (cdabs(rt) .gt. ax) tem = tem*rt 268 268 if (cdabs(rt) .gt. ax) tem = tem*rt
rt = rt+tem 269 269 rt = rt+tem
d = (h+tem)*d/h 270 270 d = (h+tem)*d/h
h = h+tem 271 271 h = h+tem
end subroutine trans_rt 272 272 end subroutine trans_rt
273 273
subroutine deflated_work(errcode) 274 274 subroutine deflated_work(errcode)
! errcode=0 => no errors 275 275 ! errcode=0 => no errors
! errcode=1 => jk>itmax or convergence of second kind achieved 276 276 ! errcode=1 => jk>itmax or convergence of second kind achieved
integer :: errcode,flag 277 277 integer :: errcode,flag
278 278
flag=1 279 279 flag=1
loop1: do while(flag==1) 280 280 loop1: do while(flag==1)
errcode=0 281 281 errcode=0
jk = jk+1 282 282 jk = jk+1
if (jk .gt. itmax) then 283 283 if (jk .gt. itmax) then
ier=33 284 284 ier=33
errcode=1 285 285 errcode=1
return 286 286 return
end if 287 287 end if
frt = f(rt) 288 288 frt = f(rt)
fprt = frt 289 289 fprt = frt
if (l /= 1) then 290 290 if (l /= 1) then
lm1 = l-1 291 291 lm1 = l-1
do i=1,lm1 292 292 do i=1,lm1
tem = rt - x(i) 293 293 tem = rt - x(i)
if (cdabs(tem) .eq. rzero) then 294 294 if (cdabs(tem) .eq. rzero) then
!if (ic .ne. 0) go to 15 !! ?? possible? 295 295 !if (ic .ne. 0) go to 15 !! ?? possible?
call trans_rt() 296 296 call trans_rt()
cycle loop1 297 297 cycle loop1