lxsd / fvn

Download as
Browse Code ยป

Commit e9d7fe24b46cc1f7f6e226744757b297f43e6828

Authored by wdaniau
1 parent f6479a913d
Exists in master and in 2 other branches geevx_2020, multi

Added fvn_lmdif to fvn_misc 

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

Showing 1 changed file with 1427 additions and 0 deletions Inline Diff

fvn_misc/fvn_misc.f90
Diff comments   View file @ e9d7fe2
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.0d0,0.0d0)/,p1/(0.1d0,0.0d0)/, & 114 114 data zero/(0.0d0,0.0d0)/,p1/(0.1d0,0.0d0)/, &
one/(1.0d0,0.0d0)/,four/(4.0d0,0.0d0)/, & 115 115 one/(1.0d0,0.0d0)/,four/(4.0d0,0.0d0)/, &
p5/(0.5d0,0.0d0)/, & 116 116 p5/(0.5d0,0.0d0)/, &
rzero/0.0d0/,rten/10.0d0/,rhun/100.0d0/, & 117 117 rzero/0.0d0/,rten/10.0d0/,rhun/100.0d0/, &
ax/0.1d0/,ickmax/3/,rp01/0.01d0/ 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 = sqrt(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
end if 298 298 end if
fprt = fprt/tem 299 299 fprt = fprt/tem
end do 300 300 end do
end if 301 301 end if
flag=0 302 302 flag=0
end do loop1 303 303 end do loop1
304 304
if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then 305 305 if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then
errcode=1 306 306 errcode=1
return 307 307 return
end if 308 308 end if
309 309
end subroutine deflated_work 310 310 end subroutine deflated_work
311 311
end subroutine 312 312 end subroutine
313 313
314 !
315 !
316 !
317 ! Non linear least square using Levenberg-Marquardt algorithm and
318 ! a finite difference jacobian
319 !
320 ! This uses MINPACK Routines (http://www.netlib.org/minpack)
321 ! Converted to fortran90 by Alan Miller amiller @ bigpond.net.au
322 !
323
324 subroutine fvn_lmdif(zef,m,n,a,info,tol)
325 integer(4), intent(in) :: m
326 integer(4), intent(in) :: n
327 real(8), dimension(:), intent(inout) :: a
328 integer(4), intent(out) :: info
329 real(8), intent(in), optional :: tol
330
331 real(8) :: rtol
332 real(8), dimension(:), allocatable :: fvec
333 integer(4), dimension(:), allocatable :: iwa
334
335 interface
336 subroutine zef(m,n,a,fvec,iflag)
337 integer(4), intent(in) :: m
338 integer(4), intent(in) :: n
339 real(8), dimension(:), intent(in) :: a
340 real(8), dimension(:), intent(inout) :: fvec
341 integer(4), intent(inout) :: iflag
342 end subroutine
343 end interface
344
345 integer(4) :: maxfev, mode, nfev, nprint
346 real(8) :: epsfcn, ftol, gtol, xtol, fjac(m,n)
347 real(8), parameter :: factor = 100._8, zero = 0.0_8
348
349 allocate(fvec(m),iwa(n))
350
351 rtol=sqrt(epsilon(0.d0))
352 if (present(tol)) rtol=tol
353
354 info = 0
355
356 ! check the input parameters for errors.
357
358 if (n <= 0 .or. m < n .or. rtol < zero) return
359
360 ! call lmdif.
361
362 maxfev = 200*(n + 1)
363 ftol = rtol
364 xtol = rtol
365 gtol = zero
366 epsfcn = zero
367 mode = 1
368 nprint = 0
369
370 call lmdif(zef, m, n, a, fvec, ftol, xtol, gtol, maxfev, epsfcn, &
371 mode, factor, nprint, info, nfev, fjac, iwa)
372
373 if (info == 8) info = 4
374
375 deallocate(fvec,iwa)
376 end subroutine
377
378
379 ! MINPACK routines which are used by both LMDIF & LMDER
380 ! 25 October 2001:
381 ! Changed INTENT of iflag in several places to IN OUT.
382 ! Changed INTENT of fvec to IN OUT in user routine FCN.
383 ! Removed arguments diag and qtv from LMDIF & LMDER.
384 ! Replaced several DO loops with array operations.
385 ! amiller @ bigpond.net.au
386
387
388 ! **********
389
390 ! subroutine lmdif1
391
392 ! The purpose of lmdif1 is to minimize the sum of the squares of m nonlinear
393 ! functions in n variables by a modification of the Levenberg-Marquardt
394 ! algorithm. This is done by using the more general least-squares
395 ! solver lmdif. The user must provide a subroutine which calculates the
396 ! functions. The jacobian is then calculated by a forward-difference
397 ! approximation.
398
399 ! the subroutine statement is
400
401 ! subroutine lmdif1(fcn, m, n, x, fvec, tol, info, iwa)
402
403 ! where
404
405 ! fcn is the name of the user-supplied subroutine which calculates
406 ! the functions. fcn must be declared in an external statement in the
407 ! user calling program, and should be written as follows.
408
409 ! subroutine fcn(m, n, x, fvec, iflag)
410 ! integer m, n, iflag
411 ! REAL (dp) x(n), fvec(m)
412 ! ----------
413 ! calculate the functions at x and return this vector in fvec.
414 ! ----------
415 ! return
416 ! end
417
418 ! the value of iflag should not be changed by fcn unless
419 ! the user wants to terminate execution of lmdif1.
420 ! In this case set iflag to a negative integer.
421
422 ! m is a positive integer input variable set to the number of functions.
423
424 ! n is a positive integer input variable set to the number of variables.
425 ! n must not exceed m.
426
427 ! x is an array of length n. On input x must contain an initial estimate
428 ! of the solution vector. On output x contains the final estimate of
429 ! the solution vector.
430
431 ! fvec is an output array of length m which contains
432 ! the functions evaluated at the output x.
433
434 ! tol is a nonnegative input variable. Termination occurs when the
435 ! algorithm estimates either that the relative error in the sum of
436 ! squares is at most tol or that the relative error between x and the
437 ! solution is at most tol.
438
439 ! info is an integer output variable. If the user has terminated execution,
440 ! info is set to the (negative) value of iflag. See description of fcn.
441 ! Otherwise, info is set as follows.
442
443 ! info = 0 improper input parameters.
444
445 ! info = 1 algorithm estimates that the relative error
446 ! in the sum of squares is at most tol.
447
448 ! info = 2 algorithm estimates that the relative error
449 ! between x and the solution is at most tol.
450
451 ! info = 3 conditions for info = 1 and info = 2 both hold.
452
453 ! info = 4 fvec is orthogonal to the columns of the
454 ! jacobian to machine precision.
455
456 ! info = 5 number of calls to fcn has reached or exceeded 200*(n+1).
457
458 ! info = 6 tol is too small. no further reduction in
459 ! the sum of squares is possible.
460
461 ! info = 7 tol is too small. No further improvement in
462 ! the approximate solution x is possible.
463
464 ! iwa is an integer work array of length n.
465
466 ! wa is a work array of length lwa.
467
468 ! lwa is a positive integer input variable not less than m*n+5*n+m.
469
470 ! subprograms called
471
472 ! user-supplied ...... fcn
473
474 ! minpack-supplied ... lmdif
475
476 ! argonne national laboratory. minpack project. march 1980.
477 ! burton s. garbow, kenneth e. hillstrom, jorge j. more
478
479 ! **********
480
481 SUBROUTINE lmdif(fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn, &
482 mode, factor, nprint, info, nfev, fjac, ipvt)
483
484 ! N.B. Arguments LDFJAC, DIAG, QTF, WA1, WA2, WA3 & WA4 have been removed.
485 INTEGER, PARAMETER :: dp = 8
486 INTEGER, INTENT(IN) :: m
487 INTEGER, INTENT(IN) :: n
488 REAL (dp), INTENT(IN OUT) :: x(:)
489 REAL (dp), INTENT(OUT) :: fvec(:)
490 REAL (dp), INTENT(IN) :: ftol
491 REAL (dp), INTENT(IN) :: xtol
492 REAL (dp), INTENT(IN OUT) :: gtol
493 INTEGER, INTENT(IN OUT) :: maxfev
494 REAL (dp), INTENT(IN OUT) :: epsfcn
495 INTEGER, INTENT(IN) :: mode
496 REAL (dp), INTENT(IN) :: factor
497 INTEGER, INTENT(IN) :: nprint
498 INTEGER, INTENT(OUT) :: info
499 INTEGER, INTENT(OUT) :: nfev
500 REAL (dp), INTENT(OUT) :: fjac(:,:) ! fjac(ldfjac,n)
501 INTEGER, INTENT(OUT) :: ipvt(:)
502
503 ! EXTERNAL fcn
504
505 INTERFACE
506 SUBROUTINE fcn(m, n, x, fvec, iflag)
507 INTEGER(4), INTENT(IN) :: m, n
508 REAL (8), INTENT(IN) :: x(:)
509 REAL (8), INTENT(IN OUT) :: fvec(:)
510 INTEGER(4), INTENT(IN OUT) :: iflag
511 END SUBROUTINE fcn
512 END INTERFACE
513
514 ! **********
515
516 ! subroutine lmdif
517
518 ! The purpose of lmdif is to minimize the sum of the squares of m nonlinear
519 ! functions in n variables by a modification of the Levenberg-Marquardt
520 ! algorithm. The user must provide a subroutine which calculates the
521 ! functions. The jacobian is then calculated by a forward-difference
522 ! approximation.
523
524 ! the subroutine statement is
525
526 ! subroutine lmdif(fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn,
527 ! diag, mode, factor, nprint, info, nfev, fjac,
528 ! ldfjac, ipvt, qtf, wa1, wa2, wa3, wa4)
529
530 ! N.B. 7 of these arguments have been removed in this version.
531
532 ! where
533
534 ! fcn is the name of the user-supplied subroutine which calculates the
535 ! functions. fcn must be declared in an external statement in the user
536 ! calling program, and should be written as follows.
537
538 ! subroutine fcn(m, n, x, fvec, iflag)
539 ! integer m, n, iflag
540 ! REAL (dp) x(:), fvec(m)
541 ! ----------
542 ! calculate the functions at x and return this vector in fvec.
543 ! ----------
544 ! return
545 ! end
546
547 ! the value of iflag should not be changed by fcn unless
548 ! the user wants to terminate execution of lmdif.
549 ! in this case set iflag to a negative integer.
550
551 ! m is a positive integer input variable set to the number of functions.
552
553 ! n is a positive integer input variable set to the number of variables.
554 ! n must not exceed m.
555
556 ! x is an array of length n. On input x must contain an initial estimate
557 ! of the solution vector. On output x contains the final estimate of the
558 ! solution vector.
559
560 ! fvec is an output array of length m which contains
561 ! the functions evaluated at the output x.
562
563 ! ftol is a nonnegative input variable. Termination occurs when both the
564 ! actual and predicted relative reductions in the sum of squares are at
565 ! most ftol. Therefore, ftol measures the relative error desired
566 ! in the sum of squares.
567
568 ! xtol is a nonnegative input variable. Termination occurs when the
569 ! relative error between two consecutive iterates is at most xtol.
570 ! Therefore, xtol measures the relative error desired in the approximate
571 ! solution.
572
573 ! gtol is a nonnegative input variable. Termination occurs when the cosine
574 ! of the angle between fvec and any column of the jacobian is at most
575 ! gtol in absolute value. Therefore, gtol measures the orthogonality
576 ! desired between the function vector and the columns of the jacobian.
577
578 ! maxfev is a positive integer input variable. Termination occurs when the
579 ! number of calls to fcn is at least maxfev by the end of an iteration.
580
581 ! epsfcn is an input variable used in determining a suitable step length
582 ! for the forward-difference approximation. This approximation assumes
583 ! that the relative errors in the functions are of the order of epsfcn.
584 ! If epsfcn is less than the machine precision, it is assumed that the
585 ! relative errors in the functions are of the order of the machine
586 ! precision.
587
588 ! diag is an array of length n. If mode = 1 (see below), diag is
589 ! internally set. If mode = 2, diag must contain positive entries that
590 ! serve as multiplicative scale factors for the variables.
591
592 ! mode is an integer input variable. If mode = 1, the variables will be
593 ! scaled internally. If mode = 2, the scaling is specified by the input
594 ! diag. other values of mode are equivalent to mode = 1.
595
596 ! factor is a positive input variable used in determining the initial step
597 ! bound. This bound is set to the product of factor and the euclidean
598 ! norm of diag*x if nonzero, or else to factor itself. In most cases
599 ! factor should lie in the interval (.1,100.). 100. is a generally
600 ! recommended value.
601
602 ! nprint is an integer input variable that enables controlled printing of
603 ! iterates if it is positive. In this case, fcn is called with iflag = 0
604 ! at the beginning of the first iteration and every nprint iterations
605 ! thereafter and immediately prior to return, with x and fvec available
606 ! for printing. If nprint is not positive, no special calls
607 ! of fcn with iflag = 0 are made.
608
609 ! info is an integer output variable. If the user has terminated
610 ! execution, info is set to the (negative) value of iflag.
611 ! See description of fcn. Otherwise, info is set as follows.
612
613 ! info = 0 improper input parameters.
614
615 ! info = 1 both actual and predicted relative reductions
616 ! in the sum of squares are at most ftol.
617
618 ! info = 2 relative error between two consecutive iterates <= xtol.
619
620 ! info = 3 conditions for info = 1 and info = 2 both hold.
621
622 ! info = 4 the cosine of the angle between fvec and any column of
623 ! the Jacobian is at most gtol in absolute value.
624
625 ! info = 5 number of calls to fcn has reached or exceeded maxfev.
626
627 ! info = 6 ftol is too small. no further reduction in
628 ! the sum of squares is possible.
629
630 ! info = 7 xtol is too small. no further improvement in
631 ! the approximate solution x is possible.
632
633 ! info = 8 gtol is too small. fvec is orthogonal to the
634 ! columns of the jacobian to machine precision.
635
636 ! nfev is an integer output variable set to the number of calls to fcn.
637
638 ! fjac is an output m by n array. the upper n by n submatrix
639 ! of fjac contains an upper triangular matrix r with
640 ! diagonal elements of nonincreasing magnitude such that
641
642 ! t t t
643 ! p *(jac *jac)*p = r *r,
644
645 ! where p is a permutation matrix and jac is the final calculated
646 ! Jacobian. Column j of p is column ipvt(j) (see below) of the
647 ! identity matrix. the lower trapezoidal part of fjac contains
648 ! information generated during the computation of r.
649
650 ! ldfjac is a positive integer input variable not less than m
651 ! which specifies the leading dimension of the array fjac.
652
653 ! ipvt is an integer output array of length n. ipvt defines a permutation
654 ! matrix p such that jac*p = q*r, where jac is the final calculated
655 ! jacobian, q is orthogonal (not stored), and r is upper triangular
656 ! with diagonal elements of nonincreasing magnitude.
657 ! Column j of p is column ipvt(j) of the identity matrix.
658
659 ! qtf is an output array of length n which contains
660 ! the first n elements of the vector (q transpose)*fvec.
661
662 ! wa1, wa2, and wa3 are work arrays of length n.
663
664 ! wa4 is a work array of length m.
665
666 ! subprograms called
667
668 ! user-supplied ...... fcn
669
670 ! minpack-supplied ... dpmpar,enorm,fdjac2,lmpar,qrfac
671
672 ! fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod
673
674 ! argonne national laboratory. minpack project. march 1980.
675 ! burton s. garbow, kenneth e. hillstrom, jorge j. more
676
677 ! **********
678 INTEGER :: i, iflag, iter, j, l
679 REAL (dp) :: actred, delta, dirder, epsmch, fnorm, fnorm1, gnorm, &
680 par, pnorm, prered, ratio, sum, temp, temp1, temp2, xnorm
681 REAL (dp) :: diag(n), qtf(n), wa1(n), wa2(n), wa3(n), wa4(m)
682 REAL (dp), PARAMETER :: one = 1.0_dp, p1 = 0.1_dp, p5 = 0.5_dp, &
683 p25 = 0.25_dp, p75 = 0.75_dp, p0001 = 0.0001_dp, &
684 zero = 0.0_dp
685
686 ! epsmch is the machine precision.
687
688 epsmch = EPSILON(zero)
689
690 info = 0
691 iflag = 0
692 nfev = 0
693
694 ! check the input parameters for errors.
695
696 IF (n <= 0 .OR. m < n .OR. ftol < zero .OR. xtol < zero .OR. gtol < zero &
697 .OR. maxfev <= 0 .OR. factor <= zero) GO TO 300
698 IF (mode /= 2) GO TO 20
699 DO j = 1, n
700 IF (diag(j) <= zero) GO TO 300
701 END DO
702
703 ! evaluate the function at the starting point and calculate its norm.
704
705 20 iflag = 1
706 CALL fcn(m, n, x, fvec, iflag)
707 nfev = 1
708 IF (iflag < 0) GO TO 300
709 fnorm = enorm(m, fvec)
710
711 ! initialize levenberg-marquardt parameter and iteration counter.
712
713 par = zero
714 iter = 1
715
716 ! beginning of the outer loop.
717
718 ! calculate the jacobian matrix.
719
720 30 iflag = 2
721 CALL fdjac2(fcn, m, n, x, fvec, fjac, iflag, epsfcn)
722 nfev = nfev + n
723 IF (iflag < 0) GO TO 300
724
725 ! If requested, call fcn to enable printing of iterates.
726
727 IF (nprint <= 0) GO TO 40
728 iflag = 0
729 IF (MOD(iter-1,nprint) == 0) CALL fcn(m, n, x, fvec, iflag)
730 IF (iflag < 0) GO TO 300
731
732 ! Compute the qr factorization of the jacobian.
733
734 40 CALL qrfac(m, n, fjac, .true., ipvt, wa1, wa2)
735
736 ! On the first iteration and if mode is 1, scale according
737 ! to the norms of the columns of the initial jacobian.
738
739 IF (iter /= 1) GO TO 80
740 IF (mode == 2) GO TO 60
741 DO j = 1, n
742 diag(j) = wa2(j)
743 IF (wa2(j) == zero) diag(j) = one
744 END DO
745
746 ! On the first iteration, calculate the norm of the scaled x
747 ! and initialize the step bound delta.
748
749 60 wa3(1:n) = diag(1:n)*x(1:n)
750 xnorm = enorm(n, wa3)
751 delta = factor*xnorm
752 IF (delta == zero) delta = factor
753
754 ! Form (q transpose)*fvec and store the first n components in qtf.
755
756 80 wa4(1:m) = fvec(1:m)
757 DO j = 1, n
758 IF (fjac(j,j) == zero) GO TO 120
759 sum = DOT_PRODUCT( fjac(j:m,j), wa4(j:m) )
760 temp = -sum/fjac(j,j)
761 DO i = j, m
762 wa4(i) = wa4(i) + fjac(i,j)*temp
763 END DO
764 120 fjac(j,j) = wa1(j)
765 qtf(j) = wa4(j)
766 END DO
767
768 ! compute the norm of the scaled gradient.
769
770 gnorm = zero
771 IF (fnorm == zero) GO TO 170
772 DO j = 1, n
773 l = ipvt(j)
774 IF (wa2(l) == zero) CYCLE
775 sum = zero
776 DO i = 1, j
777 sum = sum + fjac(i,j)*(qtf(i)/fnorm)
778 END DO
779 gnorm = MAX(gnorm, ABS(sum/wa2(l)))
780 END DO
781
782 ! test for convergence of the gradient norm.
783
784 170 IF (gnorm <= gtol) info = 4
785 IF (info /= 0) GO TO 300
786
787 ! rescale if necessary.
788
789 IF (mode == 2) GO TO 200
790 DO j = 1, n
791 diag(j) = MAX(diag(j), wa2(j))
792 END DO
793
794 ! beginning of the inner loop.
795
796 ! determine the Levenberg-Marquardt parameter.
797
798 200 CALL lmpar(n, fjac, ipvt, diag, qtf, delta, par, wa1, wa2)
799
800 ! store the direction p and x + p. calculate the norm of p.
801
802 DO j = 1, n
803 wa1(j) = -wa1(j)
804 wa2(j) = x(j) + wa1(j)
805 wa3(j) = diag(j)*wa1(j)
806 END DO
807 pnorm = enorm(n, wa3)
808
809 ! on the first iteration, adjust the initial step bound.
810
811 IF (iter == 1) delta = MIN(delta, pnorm)
812
813 ! evaluate the function at x + p and calculate its norm.
814
815 iflag = 1
816 CALL fcn(m, n, wa2, wa4, iflag)
817 nfev = nfev + 1
818 IF (iflag < 0) GO TO 300
819 fnorm1 = enorm(m, wa4)
820
821 ! compute the scaled actual reduction.
822
823 actred = -one
824 IF (p1*fnorm1 < fnorm) actred = one - (fnorm1/fnorm)**2
825
826 ! Compute the scaled predicted reduction and
827 ! the scaled directional derivative.
828
829 DO j = 1, n
830 wa3(j) = zero
831 l = ipvt(j)
832 temp = wa1(l)
833 DO i = 1, j
834 wa3(i) = wa3(i) + fjac(i,j)*temp
835 END DO
836 END DO
837 temp1 = enorm(n,wa3)/fnorm
838 temp2 = (SQRT(par)*pnorm)/fnorm
839 prered = temp1**2 + temp2**2/p5
840 dirder = -(temp1**2 + temp2**2)
841
842 ! compute the ratio of the actual to the predicted reduction.
843
844 ratio = zero
845 IF (prered /= zero) ratio = actred/prered
846
847 ! update the step bound.
848
849 IF (ratio <= p25) THEN
850 IF (actred >= zero) temp = p5
851 IF (actred < zero) temp = p5*dirder/(dirder + p5*actred)
852 IF (p1*fnorm1 >= fnorm .OR. temp < p1) temp = p1
853 delta = temp*MIN(delta,pnorm/p1)
854 par = par/temp
855 ELSE
856 IF (par /= zero .AND. ratio < p75) GO TO 260
857 delta = pnorm/p5
858 par = p5*par
859 END IF
860
861 ! test for successful iteration.
862
863 260 IF (ratio < p0001) GO TO 290
864
865 ! successful iteration. update x, fvec, and their norms.
866
867 DO j = 1, n
868 x(j) = wa2(j)
869 wa2(j) = diag(j)*x(j)
870 END DO
871 fvec(1:m) = wa4(1:m)
872 xnorm = enorm(n, wa2)
873 fnorm = fnorm1
874 iter = iter + 1
875
876 ! tests for convergence.
877
878 290 IF (ABS(actred) <= ftol .AND. prered <= ftol .AND. p5*ratio <= one) info = 1
879 IF (delta <= xtol*xnorm) info = 2
880 IF (ABS(actred) <= ftol .AND. prered <= ftol &
881 .AND. p5*ratio <= one .AND. info == 2) info = 3
882 IF (info /= 0) GO TO 300
883
884 ! tests for termination and stringent tolerances.
885
886 IF (nfev >= maxfev) info = 5
887 IF (ABS(actred) <= epsmch .AND. prered <= epsmch &
888 .AND. p5*ratio <= one) info = 6
889 IF (delta <= epsmch*xnorm) info = 7
890 IF (gnorm <= epsmch) info = 8
891 IF (info /= 0) GO TO 300
892
893 ! end of the inner loop. repeat if iteration unsuccessful.
894
895 IF (ratio < p0001) GO TO 200
896
897 ! end of the outer loop.
898
899 GO TO 30
900
901 ! termination, either normal or user imposed.
902
903 300 IF (iflag < 0) info = iflag
904 iflag = 0
905 IF (nprint > 0) CALL fcn(m, n, x, fvec, iflag)
906 RETURN
907
908 ! last card of subroutine lmdif.
909
910 END SUBROUTINE lmdif
911
912
913 ! **********
914
915
916 SUBROUTINE lmpar(n, r, ipvt, diag, qtb, delta, par, x, sdiag)
917
918 ! Code converted using TO_F90 by Alan Miller
919 ! Date: 1999-12-09 Time: 12:46:12
920
921 ! N.B. Arguments LDR, WA1 & WA2 have been removed.
922
923 INTEGER, PARAMETER :: dp = 8
924 INTEGER, INTENT(IN) :: n
925 REAL (dp), INTENT(IN OUT) :: r(:,:)
926 INTEGER, INTENT(IN) :: ipvt(:)
927 REAL (dp), INTENT(IN) :: diag(:)
928 REAL (dp), INTENT(IN) :: qtb(:)
929 REAL (dp), INTENT(IN) :: delta
930 REAL (dp), INTENT(OUT) :: par
931 REAL (dp), INTENT(OUT) :: x(:)
932 REAL (dp), INTENT(OUT) :: sdiag(:)
933
934 ! **********
935
936 ! subroutine lmpar
937
938 ! Given an m by n matrix a, an n by n nonsingular diagonal matrix d,
939 ! an m-vector b, and a positive number delta, the problem is to determine a
940 ! value for the parameter par such that if x solves the system
941
942 ! a*x = b , sqrt(par)*d*x = 0 ,
943
944 ! in the least squares sense, and dxnorm is the euclidean
945 ! norm of d*x, then either par is zero and
946
947 ! (dxnorm-delta) <= 0.1*delta ,
948
949 ! or par is positive and
950
951 ! abs(dxnorm-delta) <= 0.1*delta .
952
953 ! This subroutine completes the solution of the problem if it is provided
954 ! with the necessary information from the r factorization, with column
955 ! qpivoting, of a. That is, if a*p = q*r, where p is a permutation matrix,
956 ! q has orthogonal columns, and r is an upper triangular matrix with diagonal
957 ! elements of nonincreasing magnitude, then lmpar expects the full upper
958 ! triangle of r, the permutation matrix p, and the first n components of
959 ! (q transpose)*b.
960 ! On output lmpar also provides an upper triangular matrix s such that
961
962 ! t t t
963 ! p *(a *a + par*d*d)*p = s *s .
964
965 ! s is employed within lmpar and may be of separate interest.
966
967 ! Only a few iterations are generally needed for convergence of the algorithm.
968 ! If, however, the limit of 10 iterations is reached, then the output par
969 ! will contain the best value obtained so far.
970
971 ! the subroutine statement is
972
973 ! subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag, wa1,wa2)
974
975 ! where
976
977 ! n is a positive integer input variable set to the order of r.
978
979 ! r is an n by n array. on input the full upper triangle
980 ! must contain the full upper triangle of the matrix r.
981 ! On output the full upper triangle is unaltered, and the
982 ! strict lower triangle contains the strict upper triangle
983 ! (transposed) of the upper triangular matrix s.
984
985 ! ldr is a positive integer input variable not less than n
986 ! which specifies the leading dimension of the array r.
987
988 ! ipvt is an integer input array of length n which defines the
989 ! permutation matrix p such that a*p = q*r. column j of p
990 ! is column ipvt(j) of the identity matrix.
991
992 ! diag is an input array of length n which must contain the
993 ! diagonal elements of the matrix d.
994
995 ! qtb is an input array of length n which must contain the first
996 ! n elements of the vector (q transpose)*b.
997
998 ! delta is a positive input variable which specifies an upper
999 ! bound on the euclidean norm of d*x.
1000
1001 ! par is a nonnegative variable. on input par contains an
1002 ! initial estimate of the levenberg-marquardt parameter.
1003 ! on output par contains the final estimate.
1004
1005 ! x is an output array of length n which contains the least
1006 ! squares solution of the system a*x = b, sqrt(par)*d*x = 0,
1007 ! for the output par.
1008
1009 ! sdiag is an output array of length n which contains the
1010 ! diagonal elements of the upper triangular matrix s.
1011
1012 ! wa1 and wa2 are work arrays of length n.
1013
1014 ! subprograms called
1015
1016 ! minpack-supplied ... dpmpar,enorm,qrsolv
1017
1018 ! fortran-supplied ... ABS,MAX,MIN,SQRT
1019
1020 ! argonne national laboratory. minpack project. march 1980.
1021 ! burton s. garbow, kenneth e. hillstrom, jorge j. more
1022
1023 ! **********
1024 INTEGER :: iter, j, jm1, jp1, k, l, nsing
1025 REAL (dp) :: dxnorm, dwarf, fp, gnorm, parc, parl, paru, sum, temp
1026 REAL (dp) :: wa1(n), wa2(n)
1027 REAL (dp), PARAMETER :: p1 = 0.1_dp, p001 = 0.001_dp, zero = 0.0_dp
1028
1029 ! dwarf is the smallest positive magnitude.
1030
1031 dwarf = TINY(zero)
1032
1033 ! compute and store in x the gauss-newton direction. if the
1034 ! jacobian is rank-deficient, obtain a least squares solution.
1035
1036 nsing = n
1037 DO j = 1, n
1038 wa1(j) = qtb(j)
1039 IF (r(j,j) == zero .AND. nsing == n) nsing = j - 1
1040 IF (nsing < n) wa1(j) = zero
1041 END DO
1042
1043 DO k = 1, nsing
1044 j = nsing - k + 1
1045 wa1(j) = wa1(j)/r(j,j)
1046 temp = wa1(j)
1047 jm1 = j - 1
1048 wa1(1:jm1) = wa1(1:jm1) - r(1:jm1,j)*temp
1049 END DO
1050
1051 DO j = 1, n
1052 l = ipvt(j)
1053 x(l) = wa1(j)
1054 END DO
1055
1056 ! initialize the iteration counter.
1057 ! evaluate the function at the origin, and test
1058 ! for acceptance of the gauss-newton direction.
1059
1060 iter = 0
1061 wa2(1:n) = diag(1:n)*x(1:n)
1062 dxnorm = enorm(n, wa2)
1063 fp = dxnorm - delta
1064 IF (fp <= p1*delta) GO TO 220
1065
1066 ! if the jacobian is not rank deficient, the newton
1067 ! step provides a lower bound, parl, for the zero of
1068 ! the function. Otherwise set this bound to zero.
1069
1070 parl = zero
1071 IF (nsing < n) GO TO 120
1072 DO j = 1, n
1073 l = ipvt(j)
1074 wa1(j) = diag(l)*(wa2(l)/dxnorm)
1075 END DO
1076 DO j = 1, n
1077 sum = DOT_PRODUCT( r(1:j-1,j), wa1(1:j-1) )
1078 wa1(j) = (wa1(j) - sum)/r(j,j)
1079 END DO
1080 temp = enorm(n,wa1)
1081 parl = ((fp/delta)/temp)/temp
1082
1083 ! calculate an upper bound, paru, for the zero of the function.
1084
1085 120 DO j = 1, n
1086 sum = DOT_PRODUCT( r(1:j,j), qtb(1:j) )
1087 l = ipvt(j)
1088 wa1(j) = sum/diag(l)
1089 END DO
1090 gnorm = enorm(n,wa1)
1091 paru = gnorm/delta
1092 IF (paru == zero) paru = dwarf/MIN(delta,p1)
1093
1094 ! if the input par lies outside of the interval (parl,paru),
1095 ! set par to the closer endpoint.
1096
1097 par = MAX(par,parl)
1098 par = MIN(par,paru)
1099 IF (par == zero) par = gnorm/dxnorm
1100
1101 ! beginning of an iteration.
1102
1103 150 iter = iter + 1
1104
1105 ! evaluate the function at the current value of par.
1106
1107 IF (par == zero) par = MAX(dwarf, p001*paru)
1108 temp = SQRT(par)
1109 wa1(1:n) = temp*diag(1:n)
1110 CALL qrsolv(n, r, ipvt, wa1, qtb, x, sdiag)
1111 wa2(1:n) = diag(1:n)*x(1:n)
1112 dxnorm = enorm(n, wa2)
1113 temp = fp
1114 fp = dxnorm - delta
1115
1116 ! if the function is small enough, accept the current value
1117 ! of par. also test for the exceptional cases where parl
1118 ! is zero or the number of iterations has reached 10.
1119
1120 IF (ABS(fp) <= p1*delta .OR. parl == zero .AND. fp <= temp &
1121 .AND. temp < zero .OR. iter == 10) GO TO 220
1122
1123 ! compute the newton correction.
1124
1125 DO j = 1, n
1126 l = ipvt(j)
1127 wa1(j) = diag(l)*(wa2(l)/dxnorm)
1128 END DO
1129 DO j = 1, n
1130 wa1(j) = wa1(j)/sdiag(j)
1131 temp = wa1(j)
1132 jp1 = j + 1
1133 wa1(jp1:n) = wa1(jp1:n) - r(jp1:n,j)*temp
1134 END DO
1135 temp = enorm(n,wa1)
1136 parc = ((fp/delta)/temp)/temp
1137
1138 ! depending on the sign of the function, update parl or paru.
1139
1140 IF (fp > zero) parl = MAX(parl,par)
1141 IF (fp < zero) paru = MIN(paru,par)
1142
1143 ! compute an improved estimate for par.
1144
1145 par = MAX(parl, par+parc)
1146
1147 ! end of an iteration.
1148
1149 GO TO 150
1150
1151 ! termination.
1152
1153 220 IF (iter == 0) par = zero
1154 RETURN
1155
1156 ! last card of subroutine lmpar.
1157
1158 END SUBROUTINE lmpar
1159
1160
1161
1162 SUBROUTINE qrfac(m, n, a, pivot, ipvt, rdiag, acnorm)
1163
1164 ! Code converted using TO_F90 by Alan Miller
1165 ! Date: 1999-12-09 Time: 12:46:17
1166
1167 ! N.B. Arguments LDA, LIPVT & WA have been removed.
1168 INTEGER, PARAMETER :: dp = 8
1169 INTEGER, INTENT(IN) :: m
1170 INTEGER, INTENT(IN) :: n
1171 REAL (dp), INTENT(IN OUT) :: a(:,:)
1172 LOGICAL, INTENT(IN) :: pivot
1173 INTEGER, INTENT(OUT) :: ipvt(:)
1174 REAL (dp), INTENT(OUT) :: rdiag(:)
1175 REAL (dp), INTENT(OUT) :: acnorm(:)
1176
1177 ! **********
1178
1179 ! subroutine qrfac
1180
1181 ! This subroutine uses Householder transformations with column pivoting
1182 ! (optional) to compute a qr factorization of the m by n matrix a.
1183 ! That is, qrfac determines an orthogonal matrix q, a permutation matrix p,
1184 ! and an upper trapezoidal matrix r with diagonal elements of nonincreasing
1185 ! magnitude, such that a*p = q*r. The householder transformation for
1186 ! column k, k = 1,2,...,min(m,n), is of the form
1187
1188 ! t
1189 ! i - (1/u(k))*u*u
1190
1191 ! where u has zeros in the first k-1 positions. The form of this
1192 ! transformation and the method of pivoting first appeared in the
1193 ! corresponding linpack subroutine.
1194
1195 ! the subroutine statement is
1196
1197 ! subroutine qrfac(m, n, a, lda, pivot, ipvt, lipvt, rdiag, acnorm, wa)
1198
1199 ! N.B. 3 of these arguments have been omitted in this version.
1200
1201 ! where
1202
1203 ! m is a positive integer input variable set to the number of rows of a.
1204
1205 ! n is a positive integer input variable set to the number of columns of a.
1206
1207 ! a is an m by n array. On input a contains the matrix for
1208 ! which the qr factorization is to be computed. On output
1209 ! the strict upper trapezoidal part of a contains the strict
1210 ! upper trapezoidal part of r, and the lower trapezoidal
1211 ! part of a contains a factored form of q (the non-trivial
1212 ! elements of the u vectors described above).
1213
1214 ! lda is a positive integer input variable not less than m
1215 ! which specifies the leading dimension of the array a.
1216
1217 ! pivot is a logical input variable. If pivot is set true,
1218 ! then column pivoting is enforced. If pivot is set false,
1219 ! then no column pivoting is done.
1220
1221 ! ipvt is an integer output array of length lipvt. ipvt
1222 ! defines the permutation matrix p such that a*p = q*r.
1223 ! Column j of p is column ipvt(j) of the identity matrix.
1224 ! If pivot is false, ipvt is not referenced.
1225
1226 ! lipvt is a positive integer input variable. If pivot is false,
1227 ! then lipvt may be as small as 1. If pivot is true, then
1228 ! lipvt must be at least n.
1229
1230 ! rdiag is an output array of length n which contains the
1231 ! diagonal elements of r.
1232
1233 ! acnorm is an output array of length n which contains the norms of the
1234 ! corresponding columns of the input matrix a.
1235 ! If this information is not needed, then acnorm can coincide with rdiag.
1236
1237 ! wa is a work array of length n. If pivot is false, then wa
1238 ! can coincide with rdiag.
1239
1240 ! subprograms called
1241
1242 ! minpack-supplied ... dpmpar,enorm
1243
1244 ! fortran-supplied ... MAX,SQRT,MIN
1245
1246 ! argonne national laboratory. minpack project. march 1980.
1247 ! burton s. garbow, kenneth e. hillstrom, jorge j. more
1248
1249 ! **********
1250 INTEGER :: i, j, jp1, k, kmax, minmn
1251 REAL (dp) :: ajnorm, epsmch, sum, temp, wa(n)
1252 REAL (dp), PARAMETER :: one = 1.0_dp, p05 = 0.05_dp, zero = 0.0_dp
1253
1254 ! epsmch is the machine precision.
1255
1256 epsmch = EPSILON(zero)
1257
1258 ! compute the initial column norms and initialize several arrays.
1259
1260 DO j = 1, n
1261 acnorm(j) = enorm(m,a(1:,j))
1262 rdiag(j) = acnorm(j)
1263 wa(j) = rdiag(j)
1264 IF (pivot) ipvt(j) = j
1265 END DO
1266
1267 ! Reduce a to r with Householder transformations.
1268
1269 minmn = MIN(m,n)
1270 DO j = 1, minmn
1271 IF (.NOT.pivot) GO TO 40
1272
1273 ! Bring the column of largest norm into the pivot position.
1274
1275 kmax = j
1276 DO k = j, n
1277 IF (rdiag(k) > rdiag(kmax)) kmax = k
1278 END DO
1279 IF (kmax == j) GO TO 40
1280 DO i = 1, m
1281 temp = a(i,j)
1282 a(i,j) = a(i,kmax)
1283 a(i,kmax) = temp
1284 END DO
1285 rdiag(kmax) = rdiag(j)
1286 wa(kmax) = wa(j)
1287 k = ipvt(j)
1288 ipvt(j) = ipvt(kmax)
1289 ipvt(kmax) = k
1290
1291 ! Compute the Householder transformation to reduce the
1292 ! j-th column of a to a multiple of the j-th unit vector.
1293
1294 40 ajnorm = enorm(m-j+1, a(j:,j))
1295 IF (ajnorm == zero) CYCLE
1296 IF (a(j,j) < zero) ajnorm = -ajnorm
1297 a(j:m,j) = a(j:m,j)/ajnorm
1298 a(j,j) = a(j,j) + one
1299
1300 ! Apply the transformation to the remaining columns and update the norms.
1301
1302 jp1 = j + 1
1303 DO k = jp1, n
1304 sum = DOT_PRODUCT( a(j:m,j), a(j:m,k) )
1305 temp = sum/a(j,j)
1306 a(j:m,k) = a(j:m,k) - temp*a(j:m,j)
1307 IF (.NOT.pivot .OR. rdiag(k) == zero) CYCLE
1308 temp = a(j,k)/rdiag(k)
1309 rdiag(k) = rdiag(k)*SQRT(MAX(zero, one-temp**2))
1310 IF (p05*(rdiag(k)/wa(k))**2 > epsmch) CYCLE
1311 rdiag(k) = enorm(m-j, a(jp1:,k))
1312 wa(k) = rdiag(k)
1313 END DO
1314 rdiag(j) = -ajnorm
1315 END DO
1316 RETURN
1317
1318 ! last card of subroutine qrfac.
1319
1320 END SUBROUTINE qrfac
1321
1322
1323
1324 SUBROUTINE qrsolv(n, r, ipvt, diag, qtb, x, sdiag)
1325
1326 ! N.B. Arguments LDR & WA have been removed.
1327 INTEGER, PARAMETER :: dp = 8
1328 INTEGER, INTENT(IN) :: n
1329 REAL (dp), INTENT(IN OUT) :: r(:,:)
1330 INTEGER, INTENT(IN) :: ipvt(:)
1331 REAL (dp), INTENT(IN) :: diag(:)
1332 REAL (dp), INTENT(IN) :: qtb(:)
1333 REAL (dp), INTENT(OUT) :: x(:)
1334 REAL (dp), INTENT(OUT) :: sdiag(:)
1335
1336 ! **********
1337
1338 ! subroutine qrsolv
1339
1340 ! Given an m by n matrix a, an n by n diagonal matrix d, and an m-vector b,
1341 ! the problem is to determine an x which solves the system
1342
1343 ! a*x = b , d*x = 0 ,
1344
1345 ! in the least squares sense.
1346
1347 ! This subroutine completes the solution of the problem if it is provided
1348 ! with the necessary information from the qr factorization, with column
1349 ! pivoting, of a. That is, if a*p = q*r, where p is a permutation matrix,
1350 ! q has orthogonal columns, and r is an upper triangular matrix with diagonal
1351 ! elements of nonincreasing magnitude, then qrsolv expects the full upper
1352 ! triangle of r, the permutation matrix p, and the first n components of
1353 ! (q transpose)*b. The system a*x = b, d*x = 0, is then equivalent to
1354
1355 ! t t
1356 ! r*z = q *b , p *d*p*z = 0 ,
1357
1358 ! where x = p*z. if this system does not have full rank,
1359 ! then a least squares solution is obtained. On output qrsolv
1360 ! also provides an upper triangular matrix s such that
1361
1362 ! t t t
1363 ! p *(a *a + d*d)*p = s *s .
1364
1365 ! s is computed within qrsolv and may be of separate interest.
1366
1367 ! the subroutine statement is
1368
1369 ! subroutine qrsolv(n, r, ldr, ipvt, diag, qtb, x, sdiag, wa)
1370
1371 ! N.B. Arguments LDR and WA have been removed in this version.
1372
1373 ! where
1374
1375 ! n is a positive integer input variable set to the order of r.
1376
1377 ! r is an n by n array. On input the full upper triangle must contain
1378 ! the full upper triangle of the matrix r.
1379 ! On output the full upper triangle is unaltered, and the strict lower
1380 ! triangle contains the strict upper triangle (transposed) of the
1381 ! upper triangular matrix s.
1382
1383 ! ldr is a positive integer input variable not less than n
1384 ! which specifies the leading dimension of the array r.
1385
1386 ! ipvt is an integer input array of length n which defines the
1387 ! permutation matrix p such that a*p = q*r. Column j of p
1388 ! is column ipvt(j) of the identity matrix.
1389
1390 ! diag is an input array of length n which must contain the
1391 ! diagonal elements of the matrix d.
1392
1393 ! qtb is an input array of length n which must contain the first
1394 ! n elements of the vector (q transpose)*b.
1395
1396 ! x is an output array of length n which contains the least
1397 ! squares solution of the system a*x = b, d*x = 0.
1398
1399 ! sdiag is an output array of length n which contains the
1400 ! diagonal elements of the upper triangular matrix s.
1401
1402 ! wa is a work array of length n.
1403
1404 ! subprograms called
1405
1406 ! fortran-supplied ... ABS,SQRT
1407
1408 ! argonne national laboratory. minpack project. march 1980.
1409 ! burton s. garbow, kenneth e. hillstrom, jorge j. more
1410
1411 ! **********
1412 INTEGER :: i, j, k, kp1, l, nsing
1413 REAL (dp) :: COS, cotan, qtbpj, SIN, sum, TAN, temp, wa(n)
1414 REAL (dp), PARAMETER :: p5 = 0.5_dp, p25 = 0.25_dp, zero = 0.0_dp
1415
1416 ! Copy r and (q transpose)*b to preserve input and initialize s.
1417 ! In particular, save the diagonal elements of r in x.
1418
1419 DO j = 1, n
1420 r(j:n,j) = r(j,j:n)
1421 x(j) = r(j,j)
1422 wa(j) = qtb(j)
1423 END DO
1424
1425 ! Eliminate the diagonal matrix d using a givens rotation.
1426
1427 DO j = 1, n
1428
1429 ! Prepare the row of d to be eliminated, locating the
1430 ! diagonal element using p from the qr factorization.
1431
1432 l = ipvt(j)
1433 IF (diag(l) == zero) CYCLE
1434 sdiag(j:n) = zero
1435 sdiag(j) = diag(l)
1436
1437 ! The transformations to eliminate the row of d modify only a single
1438 ! element of (q transpose)*b beyond the first n, which is initially zero.
1439
1440 qtbpj = zero
1441 DO k = j, n
1442
1443 ! Determine a givens rotation which eliminates the
1444 ! appropriate element in the current row of d.
1445
1446 IF (sdiag(k) == zero) CYCLE
1447 IF (ABS(r(k,k)) < ABS(sdiag(k))) THEN
1448 cotan = r(k,k)/sdiag(k)
1449 SIN = p5/SQRT(p25 + p25*cotan**2)
1450 COS = SIN*cotan
1451 ELSE
1452 TAN = sdiag(k)/r(k,k)
1453 COS = p5/SQRT(p25 + p25*TAN**2)
1454 SIN = COS*TAN
1455 END IF
1456
1457 ! Compute the modified diagonal element of r and
1458 ! the modified element of ((q transpose)*b,0).
1459
1460 r(k,k) = COS*r(k,k) + SIN*sdiag(k)
1461 temp = COS*wa(k) + SIN*qtbpj
1462 qtbpj = -SIN*wa(k) + COS*qtbpj
1463 wa(k) = temp
1464
1465 ! Accumulate the tranformation in the row of s.
1466
1467 kp1 = k + 1
1468 DO i = kp1, n
1469 temp = COS*r(i,k) + SIN*sdiag(i)
1470 sdiag(i) = -SIN*r(i,k) + COS*sdiag(i)
1471 r(i,k) = temp
1472 END DO
1473 END DO
1474
1475 ! Store the diagonal element of s and restore
1476 ! the corresponding diagonal element of r.
1477
1478 sdiag(j) = r(j,j)
1479 r(j,j) = x(j)
1480 END DO
1481
1482 ! Solve the triangular system for z. If the system is singular,
1483 ! then obtain a least squares solution.
1484
1485 nsing = n
1486 DO j = 1, n
1487 IF (sdiag(j) == zero .AND. nsing == n) nsing = j - 1
1488 IF (nsing < n) wa(j) = zero
1489 END DO
1490
1491 DO k = 1, nsing
1492 j = nsing - k + 1
1493 sum = DOT_PRODUCT( r(j+1:nsing,j), wa(j+1:nsing) )
1494 wa(j) = (wa(j) - sum)/sdiag(j)
1495 END DO
1496
1497 ! Permute the components of z back to components of x.
1498
1499 DO j = 1, n
1500 l = ipvt(j)
1501 x(l) = wa(j)
1502 END DO
1503 RETURN
1504
1505 ! last card of subroutine qrsolv.
1506
1507 END SUBROUTINE qrsolv
1508
1509
1510
1511 FUNCTION enorm(n,x) RESULT(fn_val)
1512
1513 ! Code converted using TO_F90 by Alan Miller
1514 ! Date: 1999-12-09 Time: 12:45:34
1515 INTEGER, PARAMETER :: dp = 8
1516 INTEGER, INTENT(IN) :: n
1517 REAL (dp), INTENT(IN) :: x(:)
1518 REAL (dp) :: fn_val
1519
1520 ! **********
1521
1522 ! function enorm
1523
1524 ! given an n-vector x, this function calculates the euclidean norm of x.
1525
1526 ! the euclidean norm is computed by accumulating the sum of squares in
1527 ! three different sums. The sums of squares for the small and large
1528 ! components are scaled so that no overflows occur. Non-destructive
1529 ! underflows are permitted. Underflows and overflows do not occur in the
1530 ! computation of the unscaled sum of squares for the intermediate
1531 ! components. The definitions of small, intermediate and large components
1532 ! depend on two constants, rdwarf and rgiant. The main restrictions on
1533 ! these constants are that rdwarf**2 not underflow and rgiant**2 not
1534 ! overflow. The constants given here are suitable for every known computer.
1535
1536 ! the function statement is
1537
1538 ! REAL (dp) function enorm(n,x)
1539
1540 ! where
1541
1542 ! n is a positive integer input variable.
1543
1544 ! x is an input array of length n.
1545
1546 ! subprograms called
1547
1548 ! fortran-supplied ... ABS,SQRT
1549
1550 ! argonne national laboratory. minpack project. march 1980.
1551 ! burton s. garbow, kenneth e. hillstrom, jorge j. more
1552
1553 ! **********
1554 INTEGER :: i
1555 REAL (dp) :: agiant, floatn, s1, s2, s3, xabs, x1max, x3max
1556 REAL (dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp, rdwarf = 3.834E-20_dp, &
1557 rgiant = 1.304E+19_dp
1558
1559 s1 = zero
1560 s2 = zero
1561 s3 = zero
1562 x1max = zero
1563 x3max = zero
1564 floatn = n
1565 agiant = rgiant/floatn
1566 DO i = 1, n
1567 xabs = ABS(x(i))
1568 IF (xabs > rdwarf .AND. xabs < agiant) GO TO 70
1569 IF (xabs <= rdwarf) GO TO 30
1570
1571 ! sum for large components.
1572
1573 IF (xabs <= x1max) GO TO 10
1574 s1 = one + s1*(x1max/xabs)**2
1575 x1max = xabs
1576 GO TO 20
1577
1578 10 s1 = s1 + (xabs/x1max)**2
1579
1580 20 GO TO 60
1581
1582 ! sum for small components.
1583
1584 30 IF (xabs <= x3max) GO TO 40
1585 s3 = one + s3*(x3max/xabs)**2
1586 x3max = xabs
1587 GO TO 60
1588
1589 40 IF (xabs /= zero) s3 = s3 + (xabs/x3max)**2
1590
1591 60 CYCLE
1592
1593 ! sum for intermediate components.
1594
1595 70 s2 = s2 + xabs**2
1596 END DO
1597
1598 ! calculation of norm.
1599
1600 IF (s1 == zero) GO TO 100
1601 fn_val = x1max*SQRT(s1 + (s2/x1max)/x1max)
1602 GO TO 120
1603
1604 100 IF (s2 == zero) GO TO 110
1605 IF (s2 >= x3max) fn_val = SQRT(s2*(one + (x3max/s2)*(x3max*s3)))
1606 IF (s2 < x3max) fn_val = SQRT(x3max*((s2/x3max) + (x3max*s3)))
1607 GO TO 120
1608
1609 110 fn_val = x3max*SQRT(s3)
1610
1611 120 RETURN
1612
1613 ! last card of function enorm.
1614
1615 END FUNCTION enorm
1616
1617
1618
1619 SUBROUTINE fdjac2(fcn, m, n, x, fvec, fjac, iflag, epsfcn)
1620
1621 ! Code converted using TO_F90 by Alan Miller
1622 ! Date: 1999-12-09 Time: 12:45:44
1623
1624 ! N.B. Arguments LDFJAC & WA have been removed.
1625 INTEGER, PARAMETER :: dp = 8
1626 INTEGER, INTENT(IN) :: m
1627 INTEGER, INTENT(IN) :: n
1628 REAL (dp), INTENT(IN OUT) :: x(n)
1629 REAL (dp), INTENT(IN) :: fvec(m)
1630 REAL (dp), INTENT(OUT) :: fjac(:,:) ! fjac(ldfjac,n)
1631 INTEGER, INTENT(IN OUT) :: iflag
1632 REAL (dp), INTENT(IN) :: epsfcn
1633
1634 INTERFACE
1635 SUBROUTINE fcn(m, n, x, fvec, iflag)
1636 INTEGER(4), INTENT(IN) :: m, n
1637 REAL (8), INTENT(IN) :: x(:)
1638 REAL (8), INTENT(IN OUT) :: fvec(:)
1639 INTEGER(4), INTENT(IN OUT) :: iflag
1640 END SUBROUTINE fcn
1641 END INTERFACE
1642
1643 ! **********
1644
1645 ! subroutine fdjac2
1646
1647 ! this subroutine computes a forward-difference approximation
1648 ! to the m by n jacobian matrix associated with a specified
1649 ! problem of m functions in n variables.
1650
1651 ! the subroutine statement is
1652
1653 ! subroutine fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa)
1654
1655 ! where
1656
1657 ! fcn is the name of the user-supplied subroutine which calculates the
1658 ! functions. fcn must be declared in an external statement in the user
1659 ! calling program, and should be written as follows.
1660
1661 ! subroutine fcn(m,n,x,fvec,iflag)
1662 ! integer m,n,iflag
1663 ! REAL (dp) x(n),fvec(m)
1664 ! ----------
1665 ! calculate the functions at x and
1666 ! return this vector in fvec.
1667 ! ----------
1668 ! return
1669 ! end
1670
1671 ! the value of iflag should not be changed by fcn unless
1672 ! the user wants to terminate execution of fdjac2.
1673 ! in this case set iflag to a negative integer.
1674
1675 ! m is a positive integer input variable set to the number of functions.
1676
1677 ! n is a positive integer input variable set to the number of variables.
1678 ! n must not exceed m.
1679
1680 ! x is an input array of length n.
1681
1682 ! fvec is an input array of length m which must contain the
1683 ! functions evaluated at x.
1684
1685 ! fjac is an output m by n array which contains the
1686 ! approximation to the jacobian matrix evaluated at x.
1687
1688 ! ldfjac is a positive integer input variable not less than m
1689 ! which specifies the leading dimension of the array fjac.
1690
1691 ! iflag is an integer variable which can be used to terminate
1692 ! the execution of fdjac2. see description of fcn.
1693
1694 ! epsfcn is an input variable used in determining a suitable step length
1695 ! for the forward-difference approximation. This approximation assumes
1696 ! that the relative errors in the functions are of the order of epsfcn.
1697 ! If epsfcn is less than the machine precision, it is assumed that the
1698 ! relative errors in the functions are of the order of the machine
1699 ! precision.
1700
1701 ! wa is a work array of length m.
1702
1703 ! subprograms called
1704
1705 ! user-supplied ...... fcn
1706
1707 ! minpack-supplied ... dpmpar
1708
1709 ! fortran-supplied ... ABS,MAX,SQRT
1710
1711 ! argonne national laboratory. minpack project. march 1980.
1712 ! burton s. garbow, kenneth e. hillstrom, jorge j. more
1713
1714 ! **********
1715 INTEGER :: j
1716 REAL (dp) :: eps, epsmch, h, temp, wa(m)
1717 REAL (dp), PARAMETER :: zero = 0.0_dp
1718