lxsd / fvn

Download as
Browse Code ยป

Commit 20ce530d4f1b4017a2f49a5d325d46fed9ece341

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

Wrote documentation in source for fvn_lm 

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

Showing 1 changed file with 122 additions and 141 deletions Inline Diff

fvn_misc/fvn_misc.f90
Diff comments   View file @ 20ce530
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
10 public
11 private :: lmdif,lmpar,qrfac,qrsolv,enorm,fdjac2
12 ! These are made private to not interfere with a
13 ! possibly linked minpack
14
contains 10 15 contains
! 11 16 !
! Muller 12 17 ! Muller
! 13 18 !
! 14 19 !
! 15 20 !
! William Daniau 2007 16 21 ! William Daniau 2007
! 17 22 !
! This routine is a fortran 90 port of Hans D. Mittelmann's routine muller.f 18 23 ! 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 24 ! http://plato.asu.edu/ftp/other_software/muller.f
! 20 25 !
! it can be used as a replacement for imsl routine dzanly with minor changes 21 26 ! it can be used as a replacement for imsl routine dzanly with minor changes
! 22 27 !
!----------------------------------------------------------------------- 23 28 !-----------------------------------------------------------------------
! 24 29 !
! purpose - zeros of an analytic complex function 25 30 ! purpose - zeros of an analytic complex function
! using the muller method with deflation 26 31 ! using the muller method with deflation
! 27 32 !
! usage - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax, 28 33 ! usage - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax,
! infer,ier) 29 34 ! infer,ier)
! 30 35 !
! arguments f - a complex function subprogram, f(z), written 31 36 ! arguments f - a complex function subprogram, f(z), written
! by the user specifying the equation whose 32 37 ! by the user specifying the equation whose
! roots are to be found. f must appear in 33 38 ! roots are to be found. f must appear in
! an external statement in the calling pro- 34 39 ! an external statement in the calling pro-
! gram. 35 40 ! gram.
! eps - 1st stopping criterion. let fp(z)=f(z)/p 36 41 ! eps - 1st stopping criterion. let fp(z)=f(z)/p
! where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1)) 37 42 ! where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1))
! and z(1),...,z(k-1) are previously found 38 43 ! and z(1),...,z(k-1) are previously found
! roots. if ((cdabs(f(z)).le.eps) .and. 39 44 ! roots. if ((cdabs(f(z)).le.eps) .and.
! (cdabs(fp(z)).le.eps)), then z is accepted 40 45 ! (cdabs(fp(z)).le.eps)), then z is accepted
! as a root. (input) 41 46 ! as a root. (input)
! eps1 - 2nd stopping criterion. a root is accepted 42 47 ! eps1 - 2nd stopping criterion. a root is accepted
! if two successive approximations to a given 43 48 ! if two successive approximations to a given
! root agree within eps1. (input) 44 49 ! root agree within eps1. (input)
! note. if either or both of the stopping 45 50 ! note. if either or both of the stopping
! criteria are fulfilled, the root is 46 51 ! criteria are fulfilled, the root is
! accepted. 47 52 ! accepted.
! kn - the number of known roots which must be stored 48 53 ! kn - the number of known roots which must be stored
! in x(1),...,x(kn), prior to entry to muller 49 54 ! in x(1),...,x(kn), prior to entry to muller
! nguess - the number of initial guesses provided. these 50 55 ! nguess - the number of initial guesses provided. these
! guesses must be stored in x(kn+1),..., 51 56 ! guesses must be stored in x(kn+1),...,
! x(kn+nguess). nguess must be set equal 52 57 ! x(kn+nguess). nguess must be set equal
! to zero if no guesses are provided. (input) 53 58 ! to zero if no guesses are provided. (input)
! n - the number of new roots to be found by 54 59 ! n - the number of new roots to be found by
! muller (input) 55 60 ! muller (input)
! x - a complex vector of length kn+n. x(1),..., 56 61 ! x - a complex vector of length kn+n. x(1),...,
! x(kn) on input must contain any known 57 62 ! x(kn) on input must contain any known
! roots. x(kn+1),..., x(kn+n) on input may, 58 63 ! roots. x(kn+1),..., x(kn+n) on input may,
! on user option, contain initial guesses for 59 64 ! on user option, contain initial guesses for
! the n new roots which are to be computed. 60 65 ! the n new roots which are to be computed.
! if the user does not provide an initial 61 66 ! if the user does not provide an initial
! guess, zero is used. 62 67 ! guess, zero is used.
! on output, x(kn+1),...,x(kn+n) contain the 63 68 ! on output, x(kn+1),...,x(kn+n) contain the
! approximate roots found by muller. 64 69 ! approximate roots found by muller.
! itmax - the maximum allowable number of iterations 65 70 ! itmax - the maximum allowable number of iterations
! per root (input) 66 71 ! per root (input)
! infer - an integer vector of length kn+n. on 67 72 ! infer - an integer vector of length kn+n. on
! output infer(j) contains the number of 68 73 ! output infer(j) contains the number of
! iterations used in finding the j-th root 69 74 ! iterations used in finding the j-th root
! when convergence was achieved. if 70 75 ! when convergence was achieved. if
! convergence was not obtained in itmax 71 76 ! convergence was not obtained in itmax
! iterations, infer(j) will be greater than 72 77 ! iterations, infer(j) will be greater than
! itmax (output). 73 78 ! itmax (output).
! ier - error parameter (output) 74 79 ! ier - error parameter (output)
! warning error 75 80 ! warning error
! ier = 33 indicates failure to converge with- 76 81 ! ier = 33 indicates failure to converge with-
! in itmax iterations for at least one of 77 82 ! in itmax iterations for at least one of
! the (n) new roots. 78 83 ! the (n) new roots.
! 79 84 !
! 80 85 !
! remarks muller always returns the last approximation for root j 81 86 ! remarks muller always returns the last approximation for root j
! in x(j). if the convergence criterion is satisfied, 82 87 ! in x(j). if the convergence criterion is satisfied,
! then infer(j) is less than or equal to itmax. if the 83 88 ! then infer(j) is less than or equal to itmax. if the
! convergence criterion is not satisified, then infer(j) 84 89 ! convergence criterion is not satisified, then infer(j)
! is set to either itmax+1 or itmax+k, with k greater 85 90 ! is set to either itmax+1 or itmax+k, with k greater
! than 1. infer(j) = itmax+1 indicates that muller did 86 91 ! than 1. infer(j) = itmax+1 indicates that muller did
! not obtain convergence in the allowed number of iter- 87 92 ! not obtain convergence in the allowed number of iter-
! ations. in this case, the user may wish to set itmax 88 93 ! ations. in this case, the user may wish to set itmax
! to a larger value. infer(j) = itmax+k means that con- 89 94 ! to a larger value. infer(j) = itmax+k means that con-
! vergence was obtained (on iteration k) for the defla- 90 95 ! vergence was obtained (on iteration k) for the defla-
! ted function 91 96 ! ted function
! fp(z) = f(z)/((z-z(1)...(z-z(j-1))) 92 97 ! fp(z) = f(z)/((z-z(1)...(z-z(j-1)))
! 93 98 !
! but failed for f(z). in this case, better initial 94 99 ! but failed for f(z). in this case, better initial
! guesses might help or, it might be necessary to relax 95 100 ! guesses might help or, it might be necessary to relax
! the convergence criterion. 96 101 ! the convergence criterion.
! 97 102 !
!----------------------------------------------------------------------- 98 103 !-----------------------------------------------------------------------
! 99 104 !
subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier) 100 105 subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)
implicit none 101 106 implicit none
double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq,eps1w 102 107 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 108 double complex :: d,dd,den,fprt,frt,h,rt,t1,t2,t3, &
tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, & 104 109 tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, &
zero,p1,one,four,p5 105 110 zero,p1,one,four,p5
106 111
double complex, external :: f 107 112 double complex, external :: f
integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, & 108 113 integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, &
knpng,jk,ick,nn,lm1,errcode 109 114 knpng,jk,ick,nn,lm1,errcode
double complex :: x(kn+n) 110 115 double complex :: x(kn+n)
integer :: infer(kn+n) 111 116 integer :: infer(kn+n)
112 117
113 118
data zero/(0.0d0,0.0d0)/,p1/(0.1d0,0.0d0)/, & 114 119 data zero/(0.0d0,0.0d0)/,p1/(0.1d0,0.0d0)/, &
one/(1.0d0,0.0d0)/,four/(4.0d0,0.0d0)/, & 115 120 one/(1.0d0,0.0d0)/,four/(4.0d0,0.0d0)/, &
p5/(0.5d0,0.0d0)/, & 116 121 p5/(0.5d0,0.0d0)/, &
rzero/0.0d0/,rten/10.0d0/,rhun/100.0d0/, & 117 122 rzero/0.0d0/,rten/10.0d0/,rhun/100.0d0/, &
ax/0.1d0/,ickmax/3/,rp01/0.01d0/ 118 123 ax/0.1d0/,ickmax/3/,rp01/0.01d0/
119 124
ier = 0 120 125 ier = 0
if (n .lt. 1) then ! What the hell are doing here then ... 121 126 if (n .lt. 1) then ! What the hell are doing here then ...
return 122 127 return
end if 123 128 end if
!eps1 = rten **(-nsig) 124 129 !eps1 = rten **(-nsig)
eps1w = min(eps1,rp01) 125 130 eps1w = min(eps1,rp01)
126 131
knp1 = kn+1 127 132 knp1 = kn+1
knpn = kn+n 128 133 knpn = kn+n
knpng = kn+nguess 129 134 knpng = kn+nguess
do i=1,knpn 130 135 do i=1,knpn
infer(i) = 0 131 136 infer(i) = 0
if (i .gt. knpng) x(i) = zero 132 137 if (i .gt. knpng) x(i) = zero
end do 133 138 end do
l= knp1 134 139 l= knp1
135 140
ic=0 136 141 ic=0
rloop: do while (l<=knpn) ! Main loop over new roots 137 142 rloop: do while (l<=knpn) ! Main loop over new roots
jk = 0 138 143 jk = 0
ick = 0 139 144 ick = 0
xl = x(l) 140 145 xl = x(l)
icloop: do 141 146 icloop: do
ic = 0 142 147 ic = 0
h = ax 143 148 h = ax
h = p1*h 144 149 h = p1*h
if (cdabs(xl) .gt. ax) h = p1*xl 145 150 if (cdabs(xl) .gt. ax) h = p1*xl
! first three points are 146 151 ! first three points are
! xl+h, xl-h, xl 147 152 ! xl+h, xl-h, xl
rt = xl+h 148 153 rt = xl+h
call deflated_work(errcode) 149 154 call deflated_work(errcode)
if (errcode == 1) then 150 155 if (errcode == 1) then
exit icloop 151 156 exit icloop
end if 152 157 end if
153 158
z0 = fprt 154 159 z0 = fprt
y0 = frt 155 160 y0 = frt
x0 = rt 156 161 x0 = rt
rt = xl-h 157 162 rt = xl-h
call deflated_work(errcode) 158 163 call deflated_work(errcode)
if (errcode == 1) then 159 164 if (errcode == 1) then
exit icloop 160 165 exit icloop
end if 161 166 end if
162 167
z1 = fprt 163 168 z1 = fprt
y1 = frt 164 169 y1 = frt
h = xl-rt 165 170 h = xl-rt
d = h/(rt-x0) 166 171 d = h/(rt-x0)
rt = xl 167 172 rt = xl
168 173
call deflated_work(errcode) 169 174 call deflated_work(errcode)
if (errcode == 1) then 170 175 if (errcode == 1) then
exit icloop 171 176 exit icloop
end if 172 177 end if
173 178
174 179
z2 = fprt 175 180 z2 = fprt
y2 = frt 176 181 y2 = frt
! begin main algorithm 177 182 ! begin main algorithm
iloop: do 178 183 iloop: do
dd = one + d 179 184 dd = one + d
t1 = z0*d*d 180 185 t1 = z0*d*d
t2 = z1*dd*dd 181 186 t2 = z1*dd*dd
xx = z2*dd 182 187 xx = z2*dd
t3 = z2*d 183 188 t3 = z2*d
bi = t1-t2+xx+t3 184 189 bi = t1-t2+xx+t3
den = bi*bi-four*(xx*t1-t3*(t2-xx)) 185 190 den = bi*bi-four*(xx*t1-t3*(t2-xx))
! use denominator of maximum amplitude 186 191 ! use denominator of maximum amplitude
t1 = sqrt(den) 187 192 t1 = sqrt(den)
qz = rhun*max(cdabs(bi),cdabs(t1)) 188 193 qz = rhun*max(cdabs(bi),cdabs(t1))
t2 = bi + t1 189 194 t2 = bi + t1
tpq = cdabs(t2)+qz 190 195 tpq = cdabs(t2)+qz
if (tpq .eq. qz) t2 = zero 191 196 if (tpq .eq. qz) t2 = zero
t3 = bi - t1 192 197 t3 = bi - t1
tpq = cdabs(t3) + qz 193 198 tpq = cdabs(t3) + qz
if (tpq .eq. qz) t3 = zero 194 199 if (tpq .eq. qz) t3 = zero
den = t2 195 200 den = t2
qz = cdabs(t3)-cdabs(t2) 196 201 qz = cdabs(t3)-cdabs(t2)
if (qz .gt. rzero) den = t3 197 202 if (qz .gt. rzero) den = t3
! test for zero denominator 198 203 ! test for zero denominator
if (cdabs(den) .eq. rzero) then 199 204 if (cdabs(den) .eq. rzero) then
call trans_rt() 200 205 call trans_rt()
call deflated_work(errcode) 201 206 call deflated_work(errcode)
if (errcode == 1) then 202 207 if (errcode == 1) then
exit icloop 203 208 exit icloop
end if 204 209 end if
z2 = fprt 205 210 z2 = fprt
y2 = frt 206 211 y2 = frt
cycle iloop 207 212 cycle iloop
end if 208 213 end if
209 214
210 215
d = -xx/den 211 216 d = -xx/den
d = d+d 212 217 d = d+d
h = d*h 213 218 h = d*h
rt = rt + h 214 219 rt = rt + h
! check convergence of the first kind 215 220 ! check convergence of the first kind
if (cdabs(h) .le. eps1w*max(cdabs(rt),ax)) then 216 221 if (cdabs(h) .le. eps1w*max(cdabs(rt),ax)) then
if (ic .ne. 0) then 217 222 if (ic .ne. 0) then
exit icloop 218 223 exit icloop
end if 219 224 end if
ic = 1 220 225 ic = 1
z0 = y1 221 226 z0 = y1
z1 = y2 222 227 z1 = y2
z2 = f(rt) 223 228 z2 = f(rt)
xl = rt 224 229 xl = rt
ick = ick+1 225 230 ick = ick+1
if (ick .le. ickmax) then 226 231 if (ick .le. ickmax) then
cycle iloop 227 232 cycle iloop
end if 228 233 end if
! warning error, itmax = maximum 229 234 ! warning error, itmax = maximum
jk = itmax + jk 230 235 jk = itmax + jk
ier = 33 231 236 ier = 33
end if 232 237 end if
if (ic .ne. 0) then 233 238 if (ic .ne. 0) then
cycle icloop 234 239 cycle icloop
end if 235 240 end if
call deflated_work(errcode) 236 241 call deflated_work(errcode)
if (errcode == 1) then 237 242 if (errcode == 1) then
exit icloop 238 243 exit icloop
end if 239 244 end if
240 245
do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero) 241 246 do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero)
! take remedial action to induce 242 247 ! take remedial action to induce
! convergence 243 248 ! convergence
d = d*p5 244 249 d = d*p5
h = h*p5 245 250 h = h*p5
rt = rt-h 246 251 rt = rt-h
call deflated_work(errcode) 247 252 call deflated_work(errcode)
if (errcode == 1) then 248 253 if (errcode == 1) then
exit icloop 249 254 exit icloop
end if 250 255 end if
end do 251 256 end do
z0 = z1 252 257 z0 = z1
z1 = z2 253 258 z1 = z2
z2 = fprt 254 259 z2 = fprt
y0 = y1 255 260 y0 = y1
y1 = y2 256 261 y1 = y2
y2 = frt 257 262 y2 = frt
end do iloop 258 263 end do iloop
end do icloop 259 264 end do icloop
x(l) = rt 260 265 x(l) = rt
infer(l) = jk 261 266 infer(l) = jk
l = l+1 262 267 l = l+1
end do rloop 263 268 end do rloop
264 269
contains 265 270 contains
subroutine trans_rt() 266 271 subroutine trans_rt()
tem = rten*eps1w 267 272 tem = rten*eps1w
if (cdabs(rt) .gt. ax) tem = tem*rt 268 273 if (cdabs(rt) .gt. ax) tem = tem*rt
rt = rt+tem 269 274 rt = rt+tem
d = (h+tem)*d/h 270 275 d = (h+tem)*d/h
h = h+tem 271 276 h = h+tem
end subroutine trans_rt 272 277 end subroutine trans_rt
273 278
subroutine deflated_work(errcode) 274 279 subroutine deflated_work(errcode)
! errcode=0 => no errors 275 280 ! errcode=0 => no errors
! errcode=1 => jk>itmax or convergence of second kind achieved 276 281 ! errcode=1 => jk>itmax or convergence of second kind achieved
integer :: errcode,flag 277 282 integer :: errcode,flag
278 283
flag=1 279 284 flag=1
loop1: do while(flag==1) 280 285 loop1: do while(flag==1)
errcode=0 281 286 errcode=0
jk = jk+1 282 287 jk = jk+1
if (jk .gt. itmax) then 283 288 if (jk .gt. itmax) then
ier=33 284 289 ier=33
errcode=1 285 290 errcode=1
return 286 291 return
end if 287 292 end if
frt = f(rt) 288 293 frt = f(rt)
fprt = frt 289 294 fprt = frt
if (l /= 1) then 290 295 if (l /= 1) then
lm1 = l-1 291 296 lm1 = l-1
do i=1,lm1 292 297 do i=1,lm1
tem = rt - x(i) 293 298 tem = rt - x(i)
if (cdabs(tem) .eq. rzero) then 294 299 if (cdabs(tem) .eq. rzero) then
!if (ic .ne. 0) go to 15 !! ?? possible? 295 300 !if (ic .ne. 0) go to 15 !! ?? possible?
call trans_rt() 296 301 call trans_rt()
cycle loop1 297 302 cycle loop1
end if 298 303 end if
fprt = fprt/tem 299 304 fprt = fprt/tem
end do 300 305 end do
end if 301 306 end if
flag=0 302 307 flag=0
end do loop1 303 308 end do loop1
304 309
if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then 305 310 if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then
errcode=1 306 311 errcode=1
return 307 312 return
end if 308 313 end if
309 314
end subroutine deflated_work 310 315 end subroutine deflated_work
311 316
end subroutine 312 317 end subroutine
313 318
! 314 319 !
! 315 320 !
! 316 321 !
! Non linear least square using Levenberg-Marquardt algorithm and 317 322 ! Non linear least square using Levenberg-Marquardt algorithm and
! a finite difference jacobian 318 323 ! a finite difference jacobian
! 319 324 !
! This uses MINPACK Routines (http://www.netlib.org/minpack) 320 325 ! This uses MINPACK Routines (http://www.netlib.org/minpack)
! Converted to fortran90 by Alan Miller amiller @ bigpond.net.au 321 326 ! Converted to fortran90 by Alan Miller amiller @ bigpond.net.au
! 322 327 !
328 ! The purpose of fvn_lm is to minimize the sum of the squares of m nonlinear
329 ! functions in n variables by a modification of the Levenberg-Marquardt
330 ! algorithm. This is done by using the more general least-squares
331 ! solver lmdif. The user must provide a subroutine which calculates the
332 ! functions. The jacobian is then calculated by a forward-difference
333 ! approximation.
334 !
335 ! call fvn_lm(fcn,m,n,a,info,tol)
336 !
337 ! fcn : fcn is the user-supplied subroutine which calculates
338 ! the functions. fcn must follow the following interface that must
339 ! be declared in the calling subroutine :
340 !
341 ! interface
342 ! subroutine fcn(m,n,a,fvec,iflag)
343 ! use fvn_common
344 ! integer(ip_kind), intent(in) :: m
345 ! integer(ip_kind), intent(in) :: n
346 ! real(dp_kind), dimension(:), intent(in) :: a
347 ! real(dp_kind), dimension(:), intent(inout) :: fvec
348 ! integer(ip_kind), intent(inout) :: iflag
349 ! end subroutine
350 ! end interface
351 !
352 ! This is the function which calculate the differences for which which square sum
353 ! will be minimized outputing this difference in vector fvec.
354 ! Parameters of fcn are as follows :
355 ! m : positive integer input variable set to the number of functions
356 ! (number of measurement points)
357 ! n : positive integer input variable set to the number of variables
358 ! (number of parameters in the function to fit)
359 ! a : vector of length n containing parameters for which fcn should
360 ! perform the calculation
361 ! fvec : vector of length m containing the resulting evaluation
362 ! iflag : integer normally not used, can be used to exit the
363 ! the algorithm by setting it to a negative value
364 !
365 ! m : positive integer input variable set to the number of functions
366 ! (number of measurement points)
367 ! n : positive integer input variable set to the number of variables
368 ! (number of parameters in the function to fit)
369 ! a : vector of length n, on input must contains an initial guess (or unity vector)
370 ! and on output the solution vector
371 ! info : is an output positive integer
372 ! info = 0 improper input parameters.
373 ! info = 1 algorithm estimates that the relative error
374 ! in the sum of squares is at most tol.
375 ! info = 2 algorithm estimates that the relative error
376 ! between x and the solution is at most tol.
377 ! info = 3 conditions for info = 1 and info = 2 both hold.
378 ! info = 4 fvec is orthogonal to the columns of the
379 ! jacobian to machine precision.
380 ! info = 5 number of calls to fcn has reached or exceeded 200*(n+1).
381 ! info = 6 tol is too small. no further reduction in
382 ! the sum of squares is possible.
383 ! info = 7 tol is too small. No further improvement in
384 ! the approximate solution x is possible.
385 ! tol : is an optional positive value. Termination occurs when the
386 ! algorithm estimates either that the relative error in the sum of
387 ! squares is at most tol or that the relative error between x and the
388 ! solution is at most tol. If not provided default value is :
389 ! sqrt(epsilon(0.0d0))
390 !
391 subroutine fvn_lm(fcn,m,n,a,info,tol)
392 integer(ip_kind), intent(in) :: m
393 integer(ip_kind), intent(in) :: n
394 real(dp_kind), dimension(:), intent(inout) :: a
395 integer(ip_kind), intent(out) :: info
396 real(dp_kind), intent(in), optional :: tol
323 397
subroutine fvn_lmdif(zef,m,n,a,info,tol) 324 398 real(dp_kind) :: rtol
integer(4), intent(in) :: m 325 399 real(dp_kind), dimension(:), allocatable :: fvec
integer(4), intent(in) :: n 326 400 integer(ip_kind), dimension(:), allocatable :: iwa
real(8), dimension(:), intent(inout) :: a 327
integer(4), intent(out) :: info 328
real(8), intent(in), optional :: tol 329
330 401
real(8) :: rtol 331
real(8), dimension(:), allocatable :: fvec 332
integer(4), dimension(:), allocatable :: iwa 333
334
interface 335 402 interface
subroutine zef(m,n,a,fvec,iflag) 336 403 subroutine fcn(m,n,a,fvec,iflag)
integer(4), intent(in) :: m 337 404 use fvn_common
integer(4), intent(in) :: n 338 405 integer(ip_kind), intent(in) :: m
real(8), dimension(:), intent(in) :: a 339 406 integer(ip_kind), intent(in) :: n
real(8), dimension(:), intent(inout) :: fvec 340 407 real(dp_kind), dimension(:), intent(in) :: a
integer(4), intent(inout) :: iflag 341 408 real(dp_kind), dimension(:), intent(inout) :: fvec
409 integer(ip_kind), intent(inout) :: iflag
end subroutine 342 410 end subroutine
end interface 343 411 end interface
344 412
integer(4) :: maxfev, mode, nfev, nprint 345 413 integer(ip_kind) :: maxfev, mode, nfev, nprint
real(8) :: epsfcn, ftol, gtol, xtol, fjac(m,n) 346 414 real(dp_kind) :: epsfcn, ftol, gtol, xtol, fjac(m,n)
real(8), parameter :: factor = 100._8, zero = 0.0_8 347 415 real(dp_kind), parameter :: factor = 100._8, zero = 0.0_8
348 416
allocate(fvec(m),iwa(n)) 349 417 allocate(fvec(m),iwa(n))
350 418
rtol=sqrt(epsilon(0.d0)) 351 419 rtol=sqrt(epsilon(0.d0))
if (present(tol)) rtol=tol 352 420 if (present(tol)) rtol=tol
353 421
info = 0 354 422 info = 0
355 423
! check the input parameters for errors. 356 424 ! check the input parameters for errors.
357 425
if (n <= 0 .or. m < n .or. rtol < zero) return 358 426 if (n <= 0 .or. m < n .or. rtol < zero) return
359 427
! call lmdif. 360 428 ! call lmdif.
361 429
maxfev = 200*(n + 1) 362 430 maxfev = 200*(n + 1)
ftol = rtol 363 431 ftol = rtol
xtol = rtol 364 432 xtol = rtol
gtol = zero 365 433 gtol = zero
epsfcn = zero 366 434 epsfcn = zero
mode = 1 367 435 mode = 1
nprint = 0 368 436 nprint = 0
369 437
call lmdif(zef, m, n, a, fvec, ftol, xtol, gtol, maxfev, epsfcn, & 370 438 call lmdif(fcn, m, n, a, fvec, ftol, xtol, gtol, maxfev, epsfcn, &
mode, factor, nprint, info, nfev, fjac, iwa) 371 439 mode, factor, nprint, info, nfev, fjac, iwa)
372 440
if (info == 8) info = 4 373 441 if (info == 8) info = 4
374 442
deallocate(fvec,iwa) 375 443 deallocate(fvec,iwa)
end subroutine 376 444 end subroutine
377 445
378 446
! MINPACK routines which are used by both LMDIF & LMDER 379 447 ! lmdif
! 25 October 2001: 380 448 ! MINPACK Subroutine (http://www.netlib.org/minpack)
! Changed INTENT of iflag in several places to IN OUT. 381 449 ! Converted to fortran90 by Alan Miller amiller @ bigpond.net.au
! Changed INTENT of fvec to IN OUT in user routine FCN. 382
! Removed arguments diag and qtv from LMDIF & LMDER. 383
! Replaced several DO loops with array operations. 384
! amiller @ bigpond.net.au 385
386
387
! ********** 388
389
! subroutine lmdif1 390
391
! The purpose of lmdif1 is to minimize the sum of the squares of m nonlinear 392
! functions in n variables by a modification of the Levenberg-Marquardt 393
! algorithm. This is done by using the more general least-squares 394
! solver lmdif. The user must provide a subroutine which calculates the 395
! functions. The jacobian is then calculated by a forward-difference 396
! approximation. 397
398
! the subroutine statement is 399
400
! subroutine lmdif1(fcn, m, n, x, fvec, tol, info, iwa) 401
402
! where 403
404
! fcn is the name of the user-supplied subroutine which calculates 405
! the functions. fcn must be declared in an external statement in the 406
! user calling program, and should be written as follows. 407
408
! subroutine fcn(m, n, x, fvec, iflag) 409
! integer m, n, iflag 410
! REAL (dp) x(n), fvec(m) 411
! ---------- 412
! calculate the functions at x and return this vector in fvec. 413
! ---------- 414
! return 415
! end 416
417
! the value of iflag should not be changed by fcn unless 418
! the user wants to terminate execution of lmdif1. 419
! In this case set iflag to a negative integer. 420
421
! m is a positive integer input variable set to the number of functions. 422
423
! n is a positive integer input variable set to the number of variables. 424
! n must not exceed m. 425
426
! x is an array of length n. On input x must contain an initial estimate 427
! of the solution vector. On output x contains the final estimate of 428
! the solution vector. 429
430
! fvec is an output array of length m which contains 431
! the functions evaluated at the output x. 432
433
! tol is a nonnegative input variable. Termination occurs when the 434
! algorithm estimates either that the relative error in the sum of 435
! squares is at most tol or that the relative error between x and the 436
! solution is at most tol. 437
438
! info is an integer output variable. If the user has terminated execution, 439
! info is set to the (negative) value of iflag. See description of fcn. 440
! Otherwise, info is set as follows. 441
442
! info = 0 improper input parameters. 443
444
! info = 1 algorithm estimates that the relative error 445
! in the sum of squares is at most tol. 446
447
! info = 2 algorithm estimates that the relative error 448
! between x and the solution is at most tol. 449
450
! info = 3 conditions for info = 1 and info = 2 both hold. 451
452
! info = 4 fvec is orthogonal to the columns of the 453
! jacobian to machine precision. 454
455
! info = 5 number of calls to fcn has reached or exceeded 200*(n+1). 456
457
! info = 6 tol is too small. no further reduction in 458
! the sum of squares is possible. 459
460
! info = 7 tol is too small. No further improvement in 461
! the approximate solution x is possible. 462
463
! iwa is an integer work array of length n. 464
465
! wa is a work array of length lwa. 466
467
! lwa is a positive integer input variable not less than m*n+5*n+m. 468
469
! subprograms called 470
471
! user-supplied ...... fcn 472
473
! minpack-supplied ... lmdif 474
475
! argonne national laboratory. minpack project. march 1980. 476
! burton s. garbow, kenneth e. hillstrom, jorge j. more 477
478
! ********** 479
480
SUBROUTINE lmdif(fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn, & 481 450 SUBROUTINE lmdif(fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn, &
mode, factor, nprint, info, nfev, fjac, ipvt) 482 451 mode, factor, nprint, info, nfev, fjac, ipvt)
483 452
! N.B. Arguments LDFJAC, DIAG, QTF, WA1, WA2, WA3 & WA4 have been removed. 484 453 ! N.B. Arguments LDFJAC, DIAG, QTF, WA1, WA2, WA3 & WA4 have been removed.
INTEGER, PARAMETER :: dp = 8 485 454 INTEGER, PARAMETER :: dp = dp_kind
INTEGER, INTENT(IN) :: m 486 455 INTEGER, INTENT(IN) :: m
INTEGER, INTENT(IN) :: n 487 456 INTEGER, INTENT(IN) :: n
REAL (dp), INTENT(IN OUT) :: x(:) 488 457 REAL (dp), INTENT(IN OUT) :: x(:)
REAL (dp), INTENT(OUT) :: fvec(:) 489 458 REAL (dp), INTENT(OUT) :: fvec(:)
REAL (dp), INTENT(IN) :: ftol 490 459 REAL (dp), INTENT(IN) :: ftol
REAL (dp), INTENT(IN) :: xtol 491 460 REAL (dp), INTENT(IN) :: xtol
REAL (dp), INTENT(IN OUT) :: gtol 492 461 REAL (dp), INTENT(IN OUT) :: gtol
INTEGER, INTENT(IN OUT) :: maxfev 493 462 INTEGER, INTENT(IN OUT) :: maxfev
REAL (dp), INTENT(IN OUT) :: epsfcn 494 463 REAL (dp), INTENT(IN OUT) :: epsfcn
INTEGER, INTENT(IN) :: mode 495 464 INTEGER, INTENT(IN) :: mode
REAL (dp), INTENT(IN) :: factor 496 465 REAL (dp), INTENT(IN) :: factor
INTEGER, INTENT(IN) :: nprint 497 466 INTEGER, INTENT(IN) :: nprint
INTEGER, INTENT(OUT) :: info 498 467 INTEGER, INTENT(OUT) :: info
INTEGER, INTENT(OUT) :: nfev 499 468 INTEGER, INTENT(OUT) :: nfev
REAL (dp), INTENT(OUT) :: fjac(:,:) ! fjac(ldfjac,n) 500 469 REAL (dp), INTENT(OUT) :: fjac(:,:) ! fjac(ldfjac,n)
INTEGER, INTENT(OUT) :: ipvt(:) 501 470 INTEGER, INTENT(OUT) :: ipvt(:)
502 471
! EXTERNAL fcn 503 472 ! EXTERNAL fcn
504 473
INTERFACE 505 474 INTERFACE
SUBROUTINE fcn(m, n, x, fvec, iflag) 506 475 SUBROUTINE fcn(m, n, x, fvec, iflag)
INTEGER(4), INTENT(IN) :: m, n 507 476 use fvn_common
REAL (8), INTENT(IN) :: x(:) 508 477 INTEGER(ip_kind), INTENT(IN) :: m, n
REAL (8), INTENT(IN OUT) :: fvec(:) 509 478 REAL (dp_kind), INTENT(IN) :: x(:)
INTEGER(4), INTENT(IN OUT) :: iflag 510 479 REAL (dp_kind), INTENT(IN OUT) :: fvec(:)
480 INTEGER(ip_kind), INTENT(IN OUT) :: iflag
END SUBROUTINE fcn 511 481 END SUBROUTINE fcn
END INTERFACE 512 482 END INTERFACE
513 483
! ********** 514 484 ! **********
515 485
! subroutine lmdif 516 486 ! subroutine lmdif
517 487
! The purpose of lmdif is to minimize the sum of the squares of m nonlinear 518 488 ! The purpose of lmdif is to minimize the sum of the squares of m nonlinear
! functions in n variables by a modification of the Levenberg-Marquardt 519 489 ! functions in n variables by a modification of the Levenberg-Marquardt
! algorithm. The user must provide a subroutine which calculates the 520 490 ! algorithm. The user must provide a subroutine which calculates the
! functions. The jacobian is then calculated by a forward-difference 521 491 ! functions. The jacobian is then calculated by a forward-difference
! approximation. 522 492 ! approximation.
523 493
! the subroutine statement is 524 494 ! the subroutine statement is
525 495
! subroutine lmdif(fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn, 526 496 ! subroutine lmdif(fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn,
! diag, mode, factor, nprint, info, nfev, fjac, 527 497 ! diag, mode, factor, nprint, info, nfev, fjac,
! ldfjac, ipvt, qtf, wa1, wa2, wa3, wa4) 528 498 ! ldfjac, ipvt, qtf, wa1, wa2, wa3, wa4)
529 499
! N.B. 7 of these arguments have been removed in this version. 530 500 ! N.B. 7 of these arguments have been removed in this version.
531 501
! where 532 502 ! where
533 503
! fcn is the name of the user-supplied subroutine which calculates the 534 504 ! fcn is the name of the user-supplied subroutine which calculates the
! functions. fcn must be declared in an external statement in the user 535 505 ! functions. fcn must be declared in an external statement in the user
! calling program, and should be written as follows. 536 506 ! calling program, and should be written as follows.
537 507
! subroutine fcn(m, n, x, fvec, iflag) 538 508 ! subroutine fcn(m, n, x, fvec, iflag)
! integer m, n, iflag 539 509 ! integer m, n, iflag
! REAL (dp) x(:), fvec(m) 540 510 ! REAL (dp) x(:), fvec(m)
! ---------- 541 511 ! ----------
! calculate the functions at x and return this vector in fvec. 542 512 ! calculate the functions at x and return this vector in fvec.
! ---------- 543 513 ! ----------
! return 544 514 ! return
! end 545 515 ! end
546 516
! the value of iflag should not be changed by fcn unless 547 517 ! the value of iflag should not be changed by fcn unless
! the user wants to terminate execution of lmdif. 548 518 ! the user wants to terminate execution of lmdif.
! in this case set iflag to a negative integer. 549 519 ! in this case set iflag to a negative integer.
550 520
! m is a positive integer input variable set to the number of functions. 551 521 ! m is a positive integer input variable set to the number of functions.
552 522
! n is a positive integer input variable set to the number of variables. 553 523 ! n is a positive integer input variable set to the number of variables.
! n must not exceed m. 554 524 ! n must not exceed m.
555 525
! x is an array of length n. On input x must contain an initial estimate 556 526 ! x is an array of length n. On input x must contain an initial estimate
! of the solution vector. On output x contains the final estimate of the 557 527 ! of the solution vector. On output x contains the final estimate of the
! solution vector. 558 528 ! solution vector.
559 529
! fvec is an output array of length m which contains 560 530 ! fvec is an output array of length m which contains
! the functions evaluated at the output x. 561 531 ! the functions evaluated at the output x.
562 532
! ftol is a nonnegative input variable. Termination occurs when both the 563 533 ! ftol is a nonnegative input variable. Termination occurs when both the
! actual and predicted relative reductions in the sum of squares are at 564 534 ! actual and predicted relative reductions in the sum of squares are at
! most ftol. Therefore, ftol measures the relative error desired 565 535 ! most ftol. Therefore, ftol measures the relative error desired
! in the sum of squares. 566 536 ! in the sum of squares.
567 537
! xtol is a nonnegative input variable. Termination occurs when the 568 538 ! xtol is a nonnegative input variable. Termination occurs when the
! relative error between two consecutive iterates is at most xtol. 569 539 ! relative error between two consecutive iterates is at most xtol.
! Therefore, xtol measures the relative error desired in the approximate 570 540 ! Therefore, xtol measures the relative error desired in the approximate
! solution. 571 541 ! solution.
572 542
! gtol is a nonnegative input variable. Termination occurs when the cosine 573 543 ! gtol is a nonnegative input variable. Termination occurs when the cosine
! of the angle between fvec and any column of the jacobian is at most 574 544 ! of the angle between fvec and any column of the jacobian is at most
! gtol in absolute value. Therefore, gtol measures the orthogonality 575 545 ! gtol in absolute value. Therefore, gtol measures the orthogonality
! desired between the function vector and the columns of the jacobian. 576 546 ! desired between the function vector and the columns of the jacobian.
577 547
! maxfev is a positive integer input variable. Termination occurs when the 578 548 ! maxfev is a positive integer input variable. Termination occurs when the
! number of calls to fcn is at least maxfev by the end of an iteration. 579 549 ! number of calls to fcn is at least maxfev by the end of an iteration.
580 550
! epsfcn is an input variable used in determining a suitable step length 581 551 ! epsfcn is an input variable used in determining a suitable step length
! for the forward-difference approximation. This approximation assumes 582 552 ! for the forward-difference approximation. This approximation assumes
! that the relative errors in the functions are of the order of epsfcn. 583 553 ! that the relative errors in the functions are of the order of epsfcn.
! If epsfcn is less than the machine precision, it is assumed that the 584 554 ! If epsfcn is less than the machine precision, it is assumed that the
! relative errors in the functions are of the order of the machine 585 555 ! relative errors in the functions are of the order of the machine
! precision. 586 556 ! precision.
587 557
! diag is an array of length n. If mode = 1 (see below), diag is 588 558 ! diag is an array of length n. If mode = 1 (see below), diag is
! internally set. If mode = 2, diag must contain positive entries that 589 559 ! internally set. If mode = 2, diag must contain positive entries that
! serve as multiplicative scale factors for the variables. 590 560 ! serve as multiplicative scale factors for the variables.
591 561
! mode is an integer input variable. If mode = 1, the variables will be 592 562 ! mode is an integer input variable. If mode = 1, the variables will be
! scaled internally. If mode = 2, the scaling is specified by the input 593 563 ! scaled internally. If mode = 2, the scaling is specified by the input
! diag. other values of mode are equivalent to mode = 1. 594 564 ! diag. other values of mode are equivalent to mode = 1.
595 565
! factor is a positive input variable used in determining the initial step 596 566 ! factor is a positive input variable used in determining the initial step
! bound. This bound is set to the product of factor and the euclidean 597 567 ! bound. This bound is set to the product of factor and the euclidean
! norm of diag*x if nonzero, or else to factor itself. In most cases 598 568 ! norm of diag*x if nonzero, or else to factor itself. In most cases
! factor should lie in the interval (.1,100.). 100. is a generally 599 569 ! factor should lie in the interval (.1,100.). 100. is a generally
! recommended value. 600 570 ! recommended value.
601 571
! nprint is an integer input variable that enables controlled printing of 602 572 ! nprint is an integer input variable that enables controlled printing of
! iterates if it is positive. In this case, fcn is called with iflag = 0 603 573 ! iterates if it is positive. In this case, fcn is called with iflag = 0
! at the beginning of the first iteration and every nprint iterations 604 574 ! at the beginning of the first iteration and every nprint iterations
! thereafter and immediately prior to return, with x and fvec available 605 575 ! thereafter and immediately prior to return, with x and fvec available
! for printing. If nprint is not positive, no special calls 606 576 ! for printing. If nprint is not positive, no special calls
! of fcn with iflag = 0 are made. 607 577 ! of fcn with iflag = 0 are made.
608 578
! info is an integer output variable. If the user has terminated 609 579 ! info is an integer output variable. If the user has terminated
! execution, info is set to the (negative) value of iflag. 610 580 ! execution, info is set to the (negative) value of iflag.
! See description of fcn. Otherwise, info is set as follows. 611 581 ! See description of fcn. Otherwise, info is set as follows.
612 582
! info = 0 improper input parameters. 613 583 ! info = 0 improper input parameters.
614 584
! info = 1 both actual and predicted relative reductions 615 585 ! info = 1 both actual and predicted relative reductions
! in the sum of squares are at most ftol. 616 586 ! in the sum of squares are at most ftol.
617 587
! info = 2 relative error between two consecutive iterates <= xtol. 618 588 ! info = 2 relative error between two consecutive iterates <= xtol.
619 589
! info = 3 conditions for info = 1 and info = 2 both hold. 620 590 ! info = 3 conditions for info = 1 and info = 2 both hold.
621 591
! info = 4 the cosine of the angle between fvec and any column of 622 592 ! info = 4 the cosine of the angle between fvec and any column of
! the Jacobian is at most gtol in absolute value. 623 593 ! the Jacobian is at most gtol in absolute value.
624 594
! info = 5 number of calls to fcn has reached or exceeded maxfev. 625 595 ! info = 5 number of calls to fcn has reached or exceeded maxfev.
626 596
! info = 6 ftol is too small. no further reduction in 627 597 ! info = 6 ftol is too small. no further reduction in
! the sum of squares is possible. 628 598 ! the sum of squares is possible.
629 599
! info = 7 xtol is too small. no further improvement in 630 600 ! info = 7 xtol is too small. no further improvement in
! the approximate solution x is possible. 631 601 ! the approximate solution x is possible.
632 602
! info = 8 gtol is too small. fvec is orthogonal to the 633 603 ! info = 8 gtol is too small. fvec is orthogonal to the
! columns of the jacobian to machine precision. 634 604 ! columns of the jacobian to machine precision.
635 605
! nfev is an integer output variable set to the number of calls to fcn. 636 606 ! nfev is an integer output variable set to the number of calls to fcn.
637 607
! fjac is an output m by n array. the upper n by n submatrix 638 608 ! fjac is an output m by n array. the upper n by n submatrix
! of fjac contains an upper triangular matrix r with 639 609 ! of fjac contains an upper triangular matrix r with
! diagonal elements of nonincreasing magnitude such that 640 610 ! diagonal elements of nonincreasing magnitude such that
641 611
! t t t 642 612 ! t t t
! p *(jac *jac)*p = r *r, 643 613 ! p *(jac *jac)*p = r *r,
644 614
! where p is a permutation matrix and jac is the final calculated 645 615 ! where p is a permutation matrix and jac is the final calculated
! Jacobian. Column j of p is column ipvt(j) (see below) of the 646 616 ! Jacobian. Column j of p is column ipvt(j) (see below) of the
! identity matrix. the lower trapezoidal part of fjac contains 647 617 ! identity matrix. the lower trapezoidal part of fjac contains
! information generated during the computation of r. 648 618 ! information generated during the computation of r.
649 619
! ldfjac is a positive integer input variable not less than m 650 620 ! ldfjac is a positive integer input variable not less than m
! which specifies the leading dimension of the array fjac. 651 621 ! which specifies the leading dimension of the array fjac.
652 622
! ipvt is an integer output array of length n. ipvt defines a permutation 653 623 ! ipvt is an integer output array of length n. ipvt defines a permutation
! matrix p such that jac*p = q*r, where jac is the final calculated 654 624 ! matrix p such that jac*p = q*r, where jac is the final calculated
! jacobian, q is orthogonal (not stored), and r is upper triangular 655 625 ! jacobian, q is orthogonal (not stored), and r is upper triangular
! with diagonal elements of nonincreasing magnitude. 656 626 ! with diagonal elements of nonincreasing magnitude.
! Column j of p is column ipvt(j) of the identity matrix. 657 627 ! Column j of p is column ipvt(j) of the identity matrix.
658 628
! qtf is an output array of length n which contains 659 629 ! qtf is an output array of length n which contains
! the first n elements of the vector (q transpose)*fvec. 660 630 ! the first n elements of the vector (q transpose)*fvec.
661 631
! wa1, wa2, and wa3 are work arrays of length n. 662 632 ! wa1, wa2, and wa3 are work arrays of length n.
663 633
! wa4 is a work array of length m. 664 634 ! wa4 is a work array of length m.
665 635
! subprograms called 666 636 ! subprograms called
667 637
! user-supplied ...... fcn 668 638 ! user-supplied ...... fcn
669 639
! minpack-supplied ... dpmpar,enorm,fdjac2,lmpar,qrfac 670 640 ! minpack-supplied ... dpmpar,enorm,fdjac2,lmpar,qrfac
671 641
! fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod 672 642 ! fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod
673 643
! argonne national laboratory. minpack project. march 1980. 674 644 ! argonne national laboratory. minpack project. march 1980.
! burton s. garbow, kenneth e. hillstrom, jorge j. more 675 645 ! burton s. garbow, kenneth e. hillstrom, jorge j. more
676 646
! ********** 677 647 ! **********
INTEGER :: i, iflag, iter, j, l 678 648 INTEGER :: i, iflag, iter, j, l
REAL (dp) :: actred, delta, dirder, epsmch, fnorm, fnorm1, gnorm, & 679 649 REAL (dp) :: actred, delta, dirder, epsmch, fnorm, fnorm1, gnorm, &
par, pnorm, prered, ratio, sum, temp, temp1, temp2, xnorm 680 650 par, pnorm, prered, ratio, sum, temp, temp1, temp2, xnorm
REAL (dp) :: diag(n), qtf(n), wa1(n), wa2(n), wa3(n), wa4(m) 681 651 REAL (dp) :: diag(n), qtf(n), wa1(n), wa2(n), wa3(n), wa4(m)
REAL (dp), PARAMETER :: one = 1.0_dp, p1 = 0.1_dp, p5 = 0.5_dp, & 682 652 REAL (dp), PARAMETER :: one = 1.0_dp, p1 = 0.1_dp, p5 = 0.5_dp, &
p25 = 0.25_dp, p75 = 0.75_dp, p0001 = 0.0001_dp, & 683 653 p25 = 0.25_dp, p75 = 0.75_dp, p0001 = 0.0001_dp, &
zero = 0.0_dp 684 654 zero = 0.0_dp
685 655
! epsmch is the machine precision. 686 656 ! epsmch is the machine precision.
687 657
epsmch = EPSILON(zero) 688 658 epsmch = EPSILON(zero)
689 659
info = 0 690 660 info = 0
iflag = 0 691 661 iflag = 0
nfev = 0 692 662 nfev = 0
693 663
! check the input parameters for errors. 694 664 ! check the input parameters for errors.
695 665
IF (n <= 0 .OR. m < n .OR. ftol < zero .OR. xtol < zero .OR. gtol < zero & 696 666 IF (n <= 0 .OR. m < n .OR. ftol < zero .OR. xtol < zero .OR. gtol < zero &
.OR. maxfev <= 0 .OR. factor <= zero) GO TO 300 697 667 .OR. maxfev <= 0 .OR. factor <= zero) GO TO 300
IF (mode /= 2) GO TO 20 698 668 IF (mode /= 2) GO TO 20
DO j = 1, n 699 669 DO j = 1, n
IF (diag(j) <= zero) GO TO 300 700 670 IF (diag(j) <= zero) GO TO 300
END DO 701 671 END DO
702 672
! evaluate the function at the starting point and calculate its norm. 703 673 ! evaluate the function at the starting point and calculate its norm.
704 674
20 iflag = 1 705 675 20 iflag = 1
CALL fcn(m, n, x, fvec, iflag) 706 676 CALL fcn(m, n, x, fvec, iflag)
nfev = 1 707 677 nfev = 1
IF (iflag < 0) GO TO 300 708 678 IF (iflag < 0) GO TO 300
fnorm = enorm(m, fvec) 709 679 fnorm = enorm(m, fvec)
710 680
! initialize levenberg-marquardt parameter and iteration counter. 711 681 ! initialize levenberg-marquardt parameter and iteration counter.
712 682
par = zero 713 683 par = zero
iter = 1 714 684 iter = 1
715 685
! beginning of the outer loop. 716 686 ! beginning of the outer loop.
717 687
! calculate the jacobian matrix. 718 688 ! calculate the jacobian matrix.
719 689
30 iflag = 2 720 690 30 iflag = 2
CALL fdjac2(fcn, m, n, x, fvec, fjac, iflag, epsfcn) 721 691 CALL fdjac2(fcn, m, n, x, fvec, fjac, iflag, epsfcn)
nfev = nfev + n 722 692 nfev = nfev + n
IF (iflag < 0) GO TO 300 723 693 IF (iflag < 0) GO TO 300
724 694
! If requested, call fcn to enable printing of iterates. 725 695 ! If requested, call fcn to enable printing of iterates.
726 696
IF (nprint <= 0) GO TO 40 727 697 IF (nprint <= 0) GO TO 40
iflag = 0 728 698 iflag = 0
IF (MOD(iter-1,nprint) == 0) CALL fcn(m, n, x, fvec, iflag) 729 699 IF (MOD(iter-1,nprint) == 0) CALL fcn(m, n, x, fvec, iflag)
IF (iflag < 0) GO TO 300 730 700 IF (iflag < 0) GO TO 300
731 701
! Compute the qr factorization of the jacobian. 732 702 ! Compute the qr factorization of the jacobian.
733 703
40 CALL qrfac(m, n, fjac, .true., ipvt, wa1, wa2) 734 704 40 CALL qrfac(m, n, fjac, .true., ipvt, wa1, wa2)
735 705
! On the first iteration and if mode is 1, scale according 736 706 ! On the first iteration and if mode is 1, scale according
! to the norms of the columns of the initial jacobian. 737 707 ! to the norms of the columns of the initial jacobian.
738 708
IF (iter /= 1) GO TO 80 739 709 IF (iter /= 1) GO TO 80
IF (mode == 2) GO TO 60 740 710 IF (mode == 2) GO TO 60
DO j = 1, n 741 711 DO j = 1, n
diag(j) = wa2(j) 742 712 diag(j) = wa2(j)
IF (wa2(j) == zero) diag(j) = one 743 713 IF (wa2(j) == zero) diag(j) = one
END DO 744 714 END DO
745 715
! On the first iteration, calculate the norm of the scaled x 746 716 ! On the first iteration, calculate the norm of the scaled x
! and initialize the step bound delta. 747 717 ! and initialize the step bound delta.
748 718
60 wa3(1:n) = diag(1:n)*x(1:n) 749 719 60 wa3(1:n) = diag(1:n)*x(1:n)
xnorm = enorm(n, wa3) 750 720 xnorm = enorm(n, wa3)
delta = factor*xnorm 751 721 delta = factor*xnorm
IF (delta == zero) delta = factor 752 722 IF (delta == zero) delta = factor
753 723
! Form (q transpose)*fvec and store the first n components in qtf. 754 724 ! Form (q transpose)*fvec and store the first n components in qtf.
755 725
80 wa4(1:m) = fvec(1:m) 756 726 80 wa4(1:m) = fvec(1:m)
DO j = 1, n 757 727 DO j = 1, n
IF (fjac(j,j) == zero) GO TO 120 758 728 IF (fjac(j,j) == zero) GO TO 120
sum = DOT_PRODUCT( fjac(j:m,j), wa4(j:m) ) 759 729 sum = DOT_PRODUCT( fjac(j:m,j), wa4(j:m) )
temp = -sum/fjac(j,j) 760 730 temp = -sum/fjac(j,j)
DO i = j, m 761 731 DO i = j, m
wa4(i) = wa4(i) + fjac(i,j)*temp 762 732 wa4(i) = wa4(i) + fjac(i,j)*temp
END DO 763 733 END DO
120 fjac(j,j) = wa1(j) 764 734 120 fjac(j,j) = wa1(j)
qtf(j) = wa4(j) 765 735 qtf(j) = wa4(j)
END DO 766 736 END DO
767 737
! compute the norm of the scaled gradient. 768 738 ! compute the norm of the scaled gradient.
769 739
gnorm = zero 770 740 gnorm = zero
IF (fnorm == zero) GO TO 170 771 741 IF (fnorm == zero) GO TO 170
DO j = 1, n 772 742 DO j = 1, n
l = ipvt(j) 773 743 l = ipvt(j)
IF (wa2(l) == zero) CYCLE 774 744 IF (wa2(l) == zero) CYCLE
sum = zero 775 745 sum = zero
DO i = 1, j 776 746 DO i = 1, j
sum = sum + fjac(i,j)*(qtf(i)/fnorm) 777 747 sum = sum + fjac(i,j)*(qtf(i)/fnorm)
END DO 778 748 END DO
gnorm = MAX(gnorm, ABS(sum/wa2(l))) 779 749 gnorm = MAX(gnorm, ABS(sum/wa2(l)))
END DO 780 750 END DO
781 751
! test for convergence of the gradient norm. 782 752 ! test for convergence of the gradient norm.
783 753
170 IF (gnorm <= gtol) info = 4 784 754 170 IF (gnorm <= gtol) info = 4
IF (info /= 0) GO TO 300 785 755 IF (info /= 0) GO TO 300
786 756
! rescale if necessary. 787 757 ! rescale if necessary.
788 758
IF (mode == 2) GO TO 200 789 759 IF (mode == 2) GO TO 200
DO j = 1, n 790 760 DO j = 1, n
diag(j) = MAX(diag(j), wa2(j)) 791 761 diag(j) = MAX(diag(j), wa2(j))
END DO 792 762 END DO
793 763
! beginning of the inner loop. 794 764 ! beginning of the inner loop.
795 765
! determine the Levenberg-Marquardt parameter. 796 766 ! determine the Levenberg-Marquardt parameter.
797 767
200 CALL lmpar(n, fjac, ipvt, diag, qtf, delta, par, wa1, wa2) 798 768 200 CALL lmpar(n, fjac, ipvt, diag, qtf, delta, par, wa1, wa2)
799 769
! store the direction p and x + p. calculate the norm of p. 800 770 ! store the direction p and x + p. calculate the norm of p.
801 771
DO j = 1, n 802 772 DO j = 1, n
wa1(j) = -wa1(j) 803 773 wa1(j) = -wa1(j)
wa2(j) = x(j) + wa1(j) 804 774 wa2(j) = x(j) + wa1(j)
wa3(j) = diag(j)*wa1(j) 805 775 wa3(j) = diag(j)*wa1(j)
END DO 806 776 END DO
pnorm = enorm(n, wa3) 807 777 pnorm = enorm(n, wa3)
808 778
! on the first iteration, adjust the initial step bound. 809 779 ! on the first iteration, adjust the initial step bound.
810 780
IF (iter == 1) delta = MIN(delta, pnorm) 811 781 IF (iter == 1) delta = MIN(delta, pnorm)
812 782
! evaluate the function at x + p and calculate its norm. 813 783 ! evaluate the function at x + p and calculate its norm.
814 784
iflag = 1 815 785 iflag = 1
CALL fcn(m, n, wa2, wa4, iflag) 816 786 CALL fcn(m, n, wa2, wa4, iflag)
nfev = nfev + 1 817 787 nfev = nfev + 1
IF (iflag < 0) GO TO 300 818 788 IF (iflag < 0) GO TO 300
fnorm1 = enorm(m, wa4) 819 789 fnorm1 = enorm(m, wa4)
820 790
! compute the scaled actual reduction. 821 791 ! compute the scaled actual reduction.
822 792
actred = -one 823 793 actred = -one
IF (p1*fnorm1 < fnorm) actred = one - (fnorm1/fnorm)**2 824 794 IF (p1*fnorm1 < fnorm) actred = one - (fnorm1/fnorm)**2
825 795
! Compute the scaled predicted reduction and 826 796 ! Compute the scaled predicted reduction and
! the scaled directional derivative. 827 797 ! the scaled directional derivative.
828 798
DO j = 1, n 829 799 DO j = 1, n
wa3(j) = zero 830 800 wa3(j) = zero
l = ipvt(j) 831 801 l = ipvt(j)
temp = wa1(l) 832 802 temp = wa1(l)
DO i = 1, j 833 803 DO i = 1, j
wa3(i) = wa3(i) + fjac(i,j)*temp 834 804 wa3(i) = wa3(i) + fjac(i,j)*temp
END DO 835 805 END DO
END DO 836 806 END DO
temp1 = enorm(n,wa3)/fnorm 837 807 temp1 = enorm(n,wa3)/fnorm
temp2 = (SQRT(par)*pnorm)/fnorm 838 808 temp2 = (SQRT(par)*pnorm)/fnorm
prered = temp1**2 + temp2**2/p5 839 809 prered = temp1**2 + temp2**2/p5
dirder = -(temp1**2 + temp2**2) 840 810 dirder = -(temp1**2 + temp2**2)
841 811
! compute the ratio of the actual to the predicted reduction. 842 812 ! compute the ratio of the actual to the predicted reduction.
843 813
ratio = zero 844 814 ratio = zero
IF (prered /= zero) ratio = actred/prered 845 815 IF (prered /= zero) ratio = actred/prered
846 816
! update the step bound. 847 817 ! update the step bound.
848 818
IF (ratio <= p25) THEN 849 819 IF (ratio <= p25) THEN
IF (actred >= zero) temp = p5 850 820 IF (actred >= zero) temp = p5
IF (actred < zero) temp = p5*dirder/(dirder + p5*actred) 851 821 IF (actred < zero) temp = p5*dirder/(dirder + p5*actred)
IF (p1*fnorm1 >= fnorm .OR. temp < p1) temp = p1 852 822 IF (p1*fnorm1 >= fnorm .OR. temp < p1) temp = p1
delta = temp*MIN(delta,pnorm/p1) 853 823 delta = temp*MIN(delta,pnorm/p1)
par = par/temp 854 824 par = par/temp
ELSE 855 825 ELSE
IF (par /= zero .AND. ratio < p75) GO TO 260 856 826 IF (par /= zero .AND. ratio < p75) GO TO 260
delta = pnorm/p5 857 827 delta = pnorm/p5
par = p5*par 858 828 par = p5*par
END IF 859 829 END IF
860 830
! test for successful iteration. 861 831 ! test for successful iteration.
862 832
260 IF (ratio < p0001) GO TO 290 863 833 260 IF (ratio < p0001) GO TO 290
864 834
! successful iteration. update x, fvec, and their norms. 865 835 ! successful iteration. update x, fvec, and their norms.
866 836
DO j = 1, n 867 837 DO j = 1, n
x(j) = wa2(j) 868 838 x(j) = wa2(j)
wa2(j) = diag(j)*x(j) 869 839 wa2(j) = diag(j)*x(j)
END DO 870 840 END DO
fvec(1:m) = wa4(1:m) 871 841 fvec(1:m) = wa4(1:m)
xnorm = enorm(n, wa2) 872 842 xnorm = enorm(n, wa2)
fnorm = fnorm1 873 843 fnorm = fnorm1
iter = iter + 1 874 844 iter = iter + 1
875 845
! tests for convergence. 876 846 ! tests for convergence.
877 847
290 IF (ABS(actred) <= ftol .AND. prered <= ftol .AND. p5*ratio <= one) info = 1 878 848 290 IF (ABS(actred) <= ftol .AND. prered <= ftol .AND. p5*ratio <= one) info = 1
IF (delta <= xtol*xnorm) info = 2 879 849 IF (delta <= xtol*xnorm) info = 2
IF (ABS(actred) <= ftol .AND. prered <= ftol & 880 850 IF (ABS(actred) <= ftol .AND. prered <= ftol &
.AND. p5*ratio <= one .AND. info == 2) info = 3 881 851 .AND. p5*ratio <= one .AND. info == 2) info = 3
IF (info /= 0) GO TO 300 882 852 IF (info /= 0) GO TO 300
883 853
! tests for termination and stringent tolerances. 884 854 ! tests for termination and stringent tolerances.
885 855
IF (nfev >= maxfev) info = 5 886 856 IF (nfev >= maxfev) info = 5
IF (ABS(actred) <= epsmch .AND. prered <= epsmch & 887 857 IF (ABS(actred) <= epsmch .AND. prered <= epsmch &
.AND. p5*ratio <= one) info = 6 888 858 .AND. p5*ratio <= one) info = 6
IF (delta <= epsmch*xnorm) info = 7 889 859 IF (delta <= epsmch*xnorm) info = 7
IF (gnorm <= epsmch) info = 8 890 860 IF (gnorm <= epsmch) info = 8
IF (info /= 0) GO TO 300 891 861 IF (info /= 0) GO TO 300
892 862
! end of the inner loop. repeat if iteration unsuccessful. 893 863 ! end of the inner loop. repeat if iteration unsuccessful.
894 864
IF (ratio < p0001) GO TO 200 895 865 IF (ratio < p0001) GO TO 200
896 866
! end of the outer loop. 897 867 ! end of the outer loop.
898 868
GO TO 30 899 869 GO TO 30
900 870
! termination, either normal or user imposed. 901 871 ! termination, either normal or user imposed.
902 872
300 IF (iflag < 0) info = iflag 903 873 300 IF (iflag < 0) info = iflag
iflag = 0 904 874 iflag = 0
IF (nprint > 0) CALL fcn(m, n, x, fvec, iflag) 905 875 IF (nprint > 0) CALL fcn(m, n, x, fvec, iflag)
RETURN 906 876 RETURN
907 877
! last card of subroutine lmdif. 908 878 ! last card of subroutine lmdif.
909 879
END SUBROUTINE lmdif 910 880 END SUBROUTINE lmdif
911 881
912 882
! ********** 913 883 ! **********
914 884
915 885 ! lmpar
886 ! MINPACK Subroutine (http://www.netlib.org/minpack)
887 ! Converted to fortran90 by Alan Miller amiller @ bigpond.net.au
SUBROUTINE lmpar(n, r, ipvt, diag, qtb, delta, par, x, sdiag) 916 888 SUBROUTINE lmpar(n, r, ipvt, diag, qtb, delta, par, x, sdiag)
917 889
! Code converted using TO_F90 by Alan Miller 918 890 ! Code converted using TO_F90 by Alan Miller
! Date: 1999-12-09 Time: 12:46:12 919 891 ! Date: 1999-12-09 Time: 12:46:12
920 892
! N.B. Arguments LDR, WA1 & WA2 have been removed. 921 893 ! N.B. Arguments LDR, WA1 & WA2 have been removed.
922 894
INTEGER, PARAMETER :: dp = 8 923 895 INTEGER, PARAMETER :: dp = dp_kind
INTEGER, INTENT(IN) :: n 924 896 INTEGER, INTENT(IN) :: n
REAL (dp), INTENT(IN OUT) :: r(:,:) 925 897 REAL (dp), INTENT(IN OUT) :: r(:,:)
INTEGER, INTENT(IN) :: ipvt(:) 926 898 INTEGER, INTENT(IN) :: ipvt(:)
REAL (dp), INTENT(IN) :: diag(:) 927 899 REAL (dp), INTENT(IN) :: diag(:)
REAL (dp), INTENT(IN) :: qtb(:) 928 900 REAL (dp), INTENT(IN) :: qtb(:)
REAL (dp), INTENT(IN) :: delta 929 901 REAL (dp), INTENT(IN) :: delta
REAL (dp), INTENT(OUT) :: par 930 902 REAL (dp), INTENT(OUT) :: par
REAL (dp), INTENT(OUT) :: x(:) 931 903 REAL (dp), INTENT(OUT) :: x(:)
REAL (dp), INTENT(OUT) :: sdiag(:) 932 904 REAL (dp), INTENT(OUT) :: sdiag(:)
933 905
! ********** 934 906 ! **********
935 907
! subroutine lmpar 936 908 ! subroutine lmpar
937 909
! Given an m by n matrix a, an n by n nonsingular diagonal matrix d, 938 910 ! Given an m by n matrix a, an n by n nonsingular diagonal matrix d,
! an m-vector b, and a positive number delta, the problem is to determine a 939 911 ! an m-vector b, and a positive number delta, the problem is to determine a
! value for the parameter par such that if x solves the system 940 912 ! value for the parameter par such that if x solves the system
941 913
! a*x = b , sqrt(par)*d*x = 0 , 942 914 ! a*x = b , sqrt(par)*d*x = 0 ,
943 915
! in the least squares sense, and dxnorm is the euclidean 944 916 ! in the least squares sense, and dxnorm is the euclidean
! norm of d*x, then either par is zero and 945 917 ! norm of d*x, then either par is zero and
946 918
! (dxnorm-delta) <= 0.1*delta , 947 919 ! (dxnorm-delta) <= 0.1*delta ,
948 920
! or par is positive and 949 921 ! or par is positive and
950 922
! abs(dxnorm-delta) <= 0.1*delta . 951 923 ! abs(dxnorm-delta) <= 0.1*delta .
952 924
! This subroutine completes the solution of the problem if it is provided 953 925 ! This subroutine completes the solution of the problem if it is provided
! with the necessary information from the r factorization, with column 954 926 ! with the necessary information from the r factorization, with column
! qpivoting, of a. That is, if a*p = q*r, where p is a permutation matrix, 955 927 ! qpivoting, of a. That is, if a*p = q*r, where p is a permutation matrix,
! q has orthogonal columns, and r is an upper triangular matrix with diagonal 956 928 ! q has orthogonal columns, and r is an upper triangular matrix with diagonal
! elements of nonincreasing magnitude, then lmpar expects the full upper 957 929 ! elements of nonincreasing magnitude, then lmpar expects the full upper
! triangle of r, the permutation matrix p, and the first n components of 958 930 ! triangle of r, the permutation matrix p, and the first n components of
! (q transpose)*b. 959 931 ! (q transpose)*b.
! On output lmpar also provides an upper triangular matrix s such that 960 932 ! On output lmpar also provides an upper triangular matrix s such that
961 933
! t t t 962 934 ! t t t
! p *(a *a + par*d*d)*p = s *s . 963 935 ! p *(a *a + par*d*d)*p = s *s .
964 936
! s is employed within lmpar and may be of separate interest. 965 937 ! s is employed within lmpar and may be of separate interest.
966 938
! Only a few iterations are generally needed for convergence of the algorithm. 967 939 ! Only a few iterations are generally needed for convergence of the algorithm.
! If, however, the limit of 10 iterations is reached, then the output par 968 940 ! If, however, the limit of 10 iterations is reached, then the output par
! will contain the best value obtained so far. 969 941 ! will contain the best value obtained so far.
970 942
! the subroutine statement is 971 943 ! the subroutine statement is
972 944
! subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag, wa1,wa2) 973 945 ! subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag, wa1,wa2)
974 946
! where 975 947 ! where
976 948
! n is a positive integer input variable set to the order of r. 977 949 ! n is a positive integer input variable set to the order of r.
978 950
! r is an n by n array. on input the full upper triangle 979 951 ! r is an n by n array. on input the full upper triangle
! must contain the full upper triangle of the matrix r. 980 952 ! must contain the full upper triangle of the matrix r.
! On output the full upper triangle is unaltered, and the 981 953 ! On output the full upper triangle is unaltered, and the
! strict lower triangle contains the strict upper triangle 982 954 ! strict lower triangle contains the strict upper triangle
! (transposed) of the upper triangular matrix s. 983 955 ! (transposed) of the upper triangular matrix s.
984 956
! ldr is a positive integer input variable not less than n 985 957 ! ldr is a positive integer input variable not less than n
! which specifies the leading dimension of the array r. 986 958 ! which specifies the leading dimension of the array r.
987 959
! ipvt is an integer input array of length n which defines the 988 960 ! ipvt is an integer input array of length n which defines the
! permutation matrix p such that a*p = q*r. column j of p 989 961 ! permutation matrix p such that a*p = q*r. column j of p
! is column ipvt(j) of the identity matrix. 990 962 ! is column ipvt(j) of the identity matrix.
991 963
! diag is an input array of length n which must contain the 992 964 ! diag is an input array of length n which must contain the
! diagonal elements of the matrix d. 993 965 ! diagonal elements of the matrix d.
994 966
! qtb is an input array of length n which must contain the first 995 967 ! qtb is an input array of length n which must contain the first
! n elements of the vector (q transpose)*b. 996 968 ! n elements of the vector (q transpose)*b.
997 969
! delta is a positive input variable which specifies an upper 998 970 ! delta is a positive input variable which specifies an upper
! bound on the euclidean norm of d*x. 999 971 ! bound on the euclidean norm of d*x.
1000 972
! par is a nonnegative variable. on input par contains an 1001 973 ! par is a nonnegative variable. on input par contains an
! initial estimate of the levenberg-marquardt parameter. 1002 974 ! initial estimate of the levenberg-marquardt parameter.
! on output par contains the final estimate. 1003 975 ! on output par contains the final estimate.
1004 976
! x is an output array of length n which contains the least 1005 977 ! x is an output array of length n which contains the least
! squares solution of the system a*x = b, sqrt(par)*d*x = 0, 1006 978 ! squares solution of the system a*x = b, sqrt(par)*d*x = 0,
! for the output par. 1007 979 ! for the output par.
1008 980
! sdiag is an output array of length n which contains the 1009 981 ! sdiag is an output array of length n which contains the
! diagonal elements of the upper triangular matrix s. 1010 982 ! diagonal elements of the upper triangular matrix s.
1011 983
! wa1 and wa2 are work arrays of length n. 1012 984 ! wa1 and wa2 are work arrays of length n.
1013 985
! subprograms called 1014 986 ! subprograms called
1015 987
! minpack-supplied ... dpmpar,enorm,qrsolv 1016 988 ! minpack-supplied ... dpmpar,enorm,qrsolv
1017 989
! fortran-supplied ... ABS,MAX,MIN,SQRT 1018 990 ! fortran-supplied ... ABS,MAX,MIN,SQRT
1019 991
! argonne national laboratory. minpack project. march 1980. 1020 992 ! argonne national laboratory. minpack project. march 1980.
! burton s. garbow, kenneth e. hillstrom, jorge j. more 1021 993 ! burton s. garbow, kenneth e. hillstrom, jorge j. more
1022 994
! ********** 1023 995 ! **********
INTEGER :: iter, j, jm1, jp1, k, l, nsing 1024 996 INTEGER :: iter, j, jm1, jp1, k, l, nsing
REAL (dp) :: dxnorm, dwarf, fp, gnorm, parc, parl, paru, sum, temp 1025 997 REAL (dp) :: dxnorm, dwarf, fp, gnorm, parc, parl, paru, sum, temp
REAL (dp) :: wa1(n), wa2(n) 1026 998 REAL (dp) :: wa1(n), wa2(n)
REAL (dp), PARAMETER :: p1 = 0.1_dp, p001 = 0.001_dp, zero = 0.0_dp 1027 999 REAL (dp), PARAMETER :: p1 = 0.1_dp, p001 = 0.001_dp, zero = 0.0_dp
1028 1000
! dwarf is the smallest positive magnitude. 1029 1001 ! dwarf is the smallest positive magnitude.
1030 1002
dwarf = TINY(zero) 1031 1003 dwarf = TINY(zero)
1032 1004
! compute and store in x the gauss-newton direction. if the 1033 1005 ! compute and store in x the gauss-newton direction. if the
! jacobian is rank-deficient, obtain a least squares solution. 1034 1006 ! jacobian is rank-deficient, obtain a least squares solution.
1035 1007
nsing = n 1036 1008 nsing = n
DO j = 1, n 1037 1009 DO j = 1, n
wa1(j) = qtb(j) 1038 1010 wa1(j) = qtb(j)
IF (r(j,j) == zero .AND. nsing == n) nsing = j - 1 1039 1011 IF (r(j,j) == zero .AND. nsing == n) nsing = j - 1
IF (nsing < n) wa1(j) = zero 1040 1012 IF (nsing < n) wa1(j) = zero
END DO 1041 1013 END DO
1042 1014
DO k = 1, nsing 1043 1015 DO k = 1, nsing
j = nsing - k + 1 1044 1016 j = nsing - k + 1
wa1(j) = wa1(j)/r(j,j) 1045 1017 wa1(j) = wa1(j)/r(j,j)
temp = wa1(j) 1046 1018 temp = wa1(j)
jm1 = j - 1 1047 1019 jm1 = j - 1
wa1(1:jm1) = wa1(1:jm1) - r(1:jm1,j)*temp 1048 1020 wa1(1:jm1) = wa1(1:jm1) - r(1:jm1,j)*temp
END DO 1049 1021 END DO
1050 1022
DO j = 1, n 1051 1023 DO j = 1, n
l = ipvt(j) 1052 1024 l = ipvt(j)
x(l) = wa1(j) 1053 1025 x(l) = wa1(j)
END DO 1054 1026 END DO
1055 1027
! initialize the iteration counter. 1056 1028 ! initialize the iteration counter.
! evaluate the function at the origin, and test 1057 1029 ! evaluate the function at the origin, and test
! for acceptance of the gauss-newton direction. 1058 1030 ! for acceptance of the gauss-newton direction.
1059 1031
iter = 0 1060 1032 iter = 0
wa2(1:n) = diag(1:n)*x(1:n) 1061 1033 wa2(1:n) = diag(1:n)*x(1:n)
dxnorm = enorm(n, wa2) 1062 1034 dxnorm = enorm(n, wa2)
fp = dxnorm - delta 1063 1035 fp = dxnorm - delta
IF (fp <= p1*delta) GO TO 220 1064 1036 IF (fp <= p1*delta) GO TO 220
1065 1037
! if the jacobian is not rank deficient, the newton 1066 1038 ! if the jacobian is not rank deficient, the newton
! step provides a lower bound, parl, for the zero of 1067 1039 ! step provides a lower bound, parl, for the zero of
! the function. Otherwise set this bound to zero. 1068 1040 ! the function. Otherwise set this bound to zero.
1069 1041
parl = zero 1070 1042 parl = zero
IF (nsing < n) GO TO 120 1071 1043 IF (nsing < n) GO TO 120
DO j = 1, n 1072 1044 DO j = 1, n
l = ipvt(j) 1073 1045 l = ipvt(j)
wa1(j) = diag(l)*(wa2(l)/dxnorm) 1074 1046 wa1(j) = diag(l)*(wa2(l)/dxnorm)
END DO 1075 1047 END DO
DO j = 1, n 1076 1048 DO j = 1, n
sum = DOT_PRODUCT( r(1:j-1,j), wa1(1:j-1) ) 1077 1049 sum = DOT_PRODUCT( r(1:j-1,j), wa1(1:j-1) )
wa1(j) = (wa1(j) - sum)/r(j,j) 1078 1050 wa1(j) = (wa1(j) - sum)/r(j,j)
END DO 1079 1051 END DO
temp = enorm(n,wa1) 1080 1052 temp = enorm(n,wa1)
parl = ((fp/delta)/temp)/temp 1081 1053 parl = ((fp/delta)/temp)/temp
1082 1054
! calculate an upper bound, paru, for the zero of the function. 1083 1055 ! calculate an upper bound, paru, for the zero of the function.
1084 1056
120 DO j = 1, n 1085 1057 120 DO j = 1, n
sum = DOT_PRODUCT( r(1:j,j), qtb(1:j) ) 1086 1058 sum = DOT_PRODUCT( r(1:j,j), qtb(1:j) )
l = ipvt(j) 1087 1059 l = ipvt(j)
wa1(j) = sum/diag(l) 1088 1060 wa1(j) = sum/diag(l)
END DO 1089 1061 END DO
gnorm = enorm(n,wa1) 1090 1062 gnorm = enorm(n,wa1)
paru = gnorm/delta 1091 1063 paru = gnorm/delta
IF (paru == zero) paru = dwarf/MIN(delta,p1) 1092 1064 IF (paru == zero) paru = dwarf/MIN(delta,p1)
1093 1065
! if the input par lies outside of the interval (parl,paru), 1094 1066 ! if the input par lies outside of the interval (parl,paru),
! set par to the closer endpoint. 1095 1067 ! set par to the closer endpoint.
1096 1068
par = MAX(par,parl) 1097 1069 par = MAX(par,parl)
par = MIN(par,paru) 1098 1070 par = MIN(par,paru)
IF (par == zero) par = gnorm/dxnorm 1099 1071 IF (par == zero) par = gnorm/dxnorm
1100 1072
! beginning of an iteration. 1101 1073 ! beginning of an iteration.
1102 1074
150 iter = iter + 1 1103 1075 150 iter = iter + 1
1104 1076
! evaluate the function at the current value of par. 1105 1077 ! evaluate the function at the current value of par.
1106 1078
IF (par == zero) par = MAX(dwarf, p001*paru) 1107 1079 IF (par == zero) par = MAX(dwarf, p001*paru)
temp = SQRT(par) 1108 1080 temp = SQRT(par)
wa1(1:n) = temp*diag(1:n) 1109 1081 wa1(1:n) = temp*diag(1:n)
CALL qrsolv(n, r, ipvt, wa1, qtb, x, sdiag) 1110 1082 CALL qrsolv(n, r, ipvt, wa1, qtb, x, sdiag)
wa2(1:n) = diag(1:n)*x(1:n) 1111 1083 wa2(1:n) = diag(1:n)*x(1:n)
dxnorm = enorm(n, wa2) 1112 1084 dxnorm = enorm(n, wa2)
temp = fp 1113 1085 temp = fp
fp = dxnorm - delta 1114 1086 fp = dxnorm - delta
1115 1087
! if the function is small enough, accept the current value 1116 1088 ! if the function is small enough, accept the current value
! of par. also test for the exceptional cases where parl 1117 1089 ! of par. also test for the exceptional cases where parl
! is zero or the number of iterations has reached 10. 1118 1090 ! is zero or the number of iterations has reached 10.
1119 1091
IF (ABS(fp) <= p1*delta .OR. parl == zero .AND. fp <= temp & 1120 1092 IF (ABS(fp) <= p1*delta .OR. parl == zero .AND. fp <= temp &
.AND. temp < zero .OR. iter == 10) GO TO 220 1121 1093 .AND. temp < zero .OR. iter == 10) GO TO 220
1122 1094
! compute the newton correction. 1123 1095 ! compute the newton correction.
1124 1096
DO j = 1, n 1125 1097 DO j = 1, n
l = ipvt(j) 1126 1098 l = ipvt(j)
wa1(j) = diag(l)*(wa2(l)/dxnorm) 1127 1099 wa1(j) = diag(l)*(wa2(l)/dxnorm)
END DO 1128 1100 END DO
DO j = 1, n 1129 1101 DO j = 1, n
wa1(j) = wa1(j)/sdiag(j) 1130 1102 wa1(j) = wa1(j)/sdiag(j)
temp = wa1(j) 1131 1103 temp = wa1(j)
jp1 = j + 1 1132 1104 jp1 = j + 1
wa1(jp1:n) = wa1(jp1:n) - r(jp1:n,j)*temp 1133 1105 wa1(jp1:n) = wa1(jp1:n) - r(jp1:n,j)*temp
END DO 1134 1106 END DO
temp = enorm(n,wa1) 1135 1107 temp = enorm(n,wa1)
parc = ((fp/delta)/temp)/temp 1136 1108 parc = ((fp/delta)/temp)/temp
1137 1109
! depending on the sign of the function, update parl or paru. 1138 1110 ! depending on the sign of the function, update parl or paru.
1139 1111
IF (fp > zero) parl = MAX(parl,par) 1140 1112 IF (fp > zero) parl = MAX(parl,par)
IF (fp < zero) paru = MIN(paru,par) 1141 1113 IF (fp < zero) paru = MIN(paru,par)
1142 1114
! compute an improved estimate for par. 1143 1115 ! compute an improved estimate for par.
1144 1116
par = MAX(parl, par+parc) 1145 1117 par = MAX(parl, par+parc)
1146 1118
! end of an iteration. 1147 1119 ! end of an iteration.
1148 1120
GO TO 150 1149 1121 GO TO 150
1150 1122
! termination. 1151 1123 ! termination.
1152 1124
220 IF (iter == 0) par = zero 1153 1125 220 IF (iter == 0) par = zero
RETURN 1154 1126 RETURN
1155 1127
! last card of subroutine lmpar. 1156 1128 ! last card of subroutine lmpar.
1157 1129
END SUBROUTINE lmpar 1158 1130 END SUBROUTINE lmpar
1159 1131
1160 1132
1161 1133 ! qrfac
1134 ! MINPACK Subroutine (http://www.netlib.org/minpack)
1135 ! Converted to fortran90 by Alan Miller amiller @ bigpond.net.au
SUBROUTINE qrfac(m, n, a, pivot, ipvt, rdiag, acnorm) 1162 1136 SUBROUTINE qrfac(m, n, a, pivot, ipvt, rdiag, acnorm)
1163 1137
! Code converted using TO_F90 by Alan Miller 1164 1138 ! Code converted using TO_F90 by Alan Miller
! Date: 1999-12-09 Time: 12:46:17 1165 1139 ! Date: 1999-12-09 Time: 12:46:17
1166 1140
! N.B. Arguments LDA, LIPVT & WA have been removed. 1167 1141 ! N.B. Arguments LDA, LIPVT & WA have been removed.
INTEGER, PARAMETER :: dp = 8 1168 1142 INTEGER, PARAMETER :: dp = dp_kind
INTEGER, INTENT(IN) :: m 1169 1143 INTEGER, INTENT(IN) :: m
INTEGER, INTENT(IN) :: n 1170 1144 INTEGER, INTENT(IN) :: n
REAL (dp), INTENT(IN OUT) :: a(:,:) 1171 1145 REAL (dp), INTENT(IN OUT) :: a(:,:)
LOGICAL, INTENT(IN) :: pivot 1172 1146 LOGICAL, INTENT(IN) :: pivot
INTEGER, INTENT(OUT) :: ipvt(:) 1173 1147 INTEGER, INTENT(OUT) :: ipvt(:)
REAL (dp), INTENT(OUT) :: rdiag(:) 1174 1148 REAL (dp), INTENT(OUT) :: rdiag(:)
REAL (dp), INTENT(OUT) :: acnorm(:) 1175 1149 REAL (dp), INTENT(OUT) :: acnorm(:)
1176 1150
! ********** 1177 1151 ! **********
1178 1152
! subroutine qrfac 1179 1153 ! subroutine qrfac
1180 1154
! This subroutine uses Householder transformations with column pivoting 1181 1155 ! This subroutine uses Householder transformations with column pivoting
! (optional) to compute a qr factorization of the m by n matrix a. 1182 1156 ! (optional) to compute a qr factorization of the m by n matrix a.
! That is, qrfac determines an orthogonal matrix q, a permutation matrix p, 1183 1157 ! That is, qrfac determines an orthogonal matrix q, a permutation matrix p,
! and an upper trapezoidal matrix r with diagonal elements of nonincreasing 1184 1158 ! and an upper trapezoidal matrix r with diagonal elements of nonincreasing
! magnitude, such that a*p = q*r. The householder transformation for 1185 1159 ! magnitude, such that a*p = q*r. The householder transformation for
! column k, k = 1,2,...,min(m,n), is of the form 1186 1160 ! column k, k = 1,2,...,min(m,n), is of the form
1187 1161
! t 1188 1162 ! t
! i - (1/u(k))*u*u 1189 1163 ! i - (1/u(k))*u*u
1190 1164
! where u has zeros in the first k-1 positions. The form of this 1191 1165 ! where u has zeros in the first k-1 positions. The form of this
! transformation and the method of pivoting first appeared in the 1192 1166 ! transformation and the method of pivoting first appeared in the
! corresponding linpack subroutine. 1193 1167 ! corresponding linpack subroutine.
1194 1168
! the subroutine statement is 1195 1169 ! the subroutine statement is
1196 1170
! subroutine qrfac(m, n, a, lda, pivot, ipvt, lipvt, rdiag, acnorm, wa) 1197 1171 ! subroutine qrfac(m, n, a, lda, pivot, ipvt, lipvt, rdiag, acnorm, wa)
1198 1172
! N.B. 3 of these arguments have been omitted in this version. 1199 1173 ! N.B. 3 of these arguments have been omitted in this version.
1200 1174
! where 1201 1175 ! where
1202 1176
! m is a positive integer input variable set to the number of rows of a. 1203 1177 ! m is a positive integer input variable set to the number of rows of a.
1204 1178
! n is a positive integer input variable set to the number of columns of a. 1205 1179 ! n is a positive integer input variable set to the number of columns of a.
1206 1180
! a is an m by n array. On input a contains the matrix for 1207 1181 ! a is an m by n array. On input a contains the matrix for
! which the qr factorization is to be computed. On output 1208 1182 ! which the qr factorization is to be computed. On output
! the strict upper trapezoidal part of a contains the strict 1209 1183 ! the strict upper trapezoidal part of a contains the strict
! upper trapezoidal part of r, and the lower trapezoidal 1210 1184 ! upper trapezoidal part of r, and the lower trapezoidal
! part of a contains a factored form of q (the non-trivial 1211 1185 ! part of a contains a factored form of q (the non-trivial
! elements of the u vectors described above). 1212 1186 ! elements of the u vectors described above).
1213 1187
! lda is a positive integer input variable not less than m 1214 1188 ! lda is a positive integer input variable not less than m
! which specifies the leading dimension of the array a. 1215 1189 ! which specifies the leading dimension of the array a.
1216 1190
! pivot is a logical input variable. If pivot is set true, 1217 1191 ! pivot is a logical input variable. If pivot is set true,
! then column pivoting is enforced. If pivot is set false, 1218 1192 ! then column pivoting is enforced. If pivot is set false,
! then no column pivoting is done. 1219 1193 ! then no column pivoting is done.
1220 1194
! ipvt is an integer output array of length lipvt. ipvt 1221 1195 ! ipvt is an integer output array of length lipvt. ipvt
! defines the permutation matrix p such that a*p = q*r. 1222 1196 ! defines the permutation matrix p such that a*p = q*r.
! Column j of p is column ipvt(j) of the identity matrix. 1223 1197 ! Column j of p is column ipvt(j) of the identity matrix.
! If pivot is false, ipvt is not referenced. 1224 1198 ! If pivot is false, ipvt is not referenced.
1225 1199
! lipvt is a positive integer input variable. If pivot is false, 1226 1200 ! lipvt is a positive integer input variable. If pivot is false,
! then lipvt may be as small as 1. If pivot is true, then 1227 1201 ! then lipvt may be as small as 1. If pivot is true, then
! lipvt must be at least n. 1228 1202 ! lipvt must be at least n.
1229 1203
! rdiag is an output array of length n which contains the 1230 1204 ! rdiag is an output array of length n which contains the
! diagonal elements of r. 1231 1205 ! diagonal elements of r.
1232 1206
! acnorm is an output array of length n which contains the norms of the 1233 1207 ! acnorm is an output array of length n which contains the norms of the
! corresponding columns of the input matrix a. 1234 1208 ! corresponding columns of the input matrix a.
! If this information is not needed, then acnorm can coincide with rdiag. 1235 1209 ! If this information is not needed, then acnorm can coincide with rdiag.
1236 1210
! wa is a work array of length n. If pivot is false, then wa 1237 1211 ! wa is a work array of length n. If pivot is false, then wa
! can coincide with rdiag. 1238 1212 ! can coincide with rdiag.
1239 1213
! subprograms called 1240 1214 ! subprograms called
1241 1215
! minpack-supplied ... dpmpar,enorm 1242 1216 ! minpack-supplied ... dpmpar,enorm
1243 1217
! fortran-supplied ... MAX,SQRT,MIN 1244 1218 ! fortran-supplied ... MAX,SQRT,MIN
1245 1219
! argonne national laboratory. minpack project. march 1980. 1246 1220 ! argonne national laboratory. minpack project. march 1980.
! burton s. garbow, kenneth e. hillstrom, jorge j. more 1247 1221 ! burton s. garbow, kenneth e. hillstrom, jorge j. more
1248 1222
! ********** 1249 1223 ! **********
INTEGER :: i, j, jp1, k, kmax, minmn 1250 1224 INTEGER :: i, j, jp1, k, kmax, minmn
REAL (dp) :: ajnorm, epsmch, sum, temp, wa(n) 1251 1225 REAL (dp) :: ajnorm, epsmch, sum, temp, wa(n)
REAL (dp), PARAMETER :: one = 1.0_dp, p05 = 0.05_dp, zero = 0.0_dp 1252 1226 REAL (dp), PARAMETER :: one = 1.0_dp, p05 = 0.05_dp, zero = 0.0_dp
1253 1227
! epsmch is the machine precision. 1254 1228 ! epsmch is the machine precision.
1255 1229
epsmch = EPSILON(zero) 1256 1230 epsmch = EPSILON(zero)
1257 1231
! compute the initial column norms and initialize several arrays. 1258 1232 ! compute the initial column norms and initialize several arrays.
1259 1233
DO j = 1, n 1260 1234 DO j = 1, n
acnorm(j) = enorm(m,a(1:,j)) 1261 1235 acnorm(j) = enorm(m,a(1:,j))
rdiag(j) = acnorm(j) 1262 1236 rdiag(j) = acnorm(j)
wa(j) = rdiag(j) 1263 1237 wa(j) = rdiag(j)
IF (pivot) ipvt(j) = j 1264 1238 IF (pivot) ipvt(j) = j
END DO 1265 1239 END DO
1266 1240
! Reduce a to r with Householder transformations. 1267 1241 ! Reduce a to r with Householder transformations.
1268 1242
minmn = MIN(m,n) 1269 1243 minmn = MIN(m,n)
DO j = 1, minmn 1270 1244 DO j = 1, minmn
IF (.NOT.pivot) GO TO 40 1271 1245 IF (.NOT.pivot) GO TO 40
1272 1246
! Bring the column of largest norm into the pivot position. 1273 1247 ! Bring the column of largest norm into the pivot position.
1274 1248
kmax = j 1275 1249 kmax = j
DO k = j, n 1276 1250 DO k = j, n
IF (rdiag(k) > rdiag(kmax)) kmax = k 1277 1251 IF (rdiag(k) > rdiag(kmax)) kmax = k
END DO 1278 1252 END DO
IF (kmax == j) GO TO 40 1279 1253 IF (kmax == j) GO TO 40
DO i = 1, m 1280 1254 DO i = 1, m
temp = a(i,j) 1281 1255 temp = a(i,j)
a(i,j) = a(i,kmax) 1282 1256 a(i,j) = a(i,kmax)
a(i,kmax) = temp 1283 1257 a(i,kmax) = temp
END DO 1284 1258 END DO
rdiag(kmax) = rdiag(j) 1285 1259 rdiag(kmax) = rdiag(j)
wa(kmax) = wa(j) 1286 1260 wa(kmax) = wa(j)
k = ipvt(j) 1287 1261 k = ipvt(j)
ipvt(j) = ipvt(kmax) 1288 1262 ipvt(j) = ipvt(kmax)
ipvt(kmax) = k 1289 1263 ipvt(kmax) = k
1290 1264
! Compute the Householder transformation to reduce the 1291 1265 ! Compute the Householder transformation to reduce the
! j-th column of a to a multiple of the j-th unit vector. 1292 1266 ! j-th column of a to a multiple of the j-th unit vector.
1293 1267
40 ajnorm = enorm(m-j+1, a(j:,j)) 1294 1268 40 ajnorm = enorm(m-j+1, a(j:,j))
IF (ajnorm == zero) CYCLE 1295 1269 IF (ajnorm == zero) CYCLE
IF (a(j,j) < zero) ajnorm = -ajnorm 1296 1270 IF (a(j,j) < zero) ajnorm = -ajnorm
a(j:m,j) = a(j:m,j)/ajnorm 1297 1271 a(j:m,j) = a(j:m,j)/ajnorm
a(j,j) = a(j,j) + one 1298 1272 a(j,j) = a(j,j) + one
1299 1273
! Apply the transformation to the remaining columns and update the norms. 1300 1274 ! Apply the transformation to the remaining columns and update the norms.
1301 1275
jp1 = j + 1 1302 1276 jp1 = j + 1
DO k = jp1, n 1303 1277 DO k = jp1, n
sum = DOT_PRODUCT( a(j:m,j), a(j:m,k) ) 1304 1278 sum = DOT_PRODUCT( a(j:m,j), a(j:m,k) )
temp = sum/a(j,j) 1305 1279 temp = sum/a(j,j)
a(j:m,k) = a(j:m,k) - temp*a(j:m,j) 1306 1280 a(j:m,k) = a(j:m,k) - temp*a(j:m,j)
IF (.NOT.pivot .OR. rdiag(k) == zero) CYCLE 1307 1281 IF (.NOT.pivot .OR. rdiag(k) == zero) CYCLE
temp = a(j,k)/rdiag(k) 1308 1282 temp = a(j,k)/rdiag(k)
rdiag(k) = rdiag(k)*SQRT(MAX(zero, one-temp**2)) 1309 1283 rdiag(k) = rdiag(k)*SQRT(MAX(zero, one-temp**2))
IF (p05*(rdiag(k)/wa(k))**2 > epsmch) CYCLE 1310 1284 IF (p05*(rdiag(k)/wa(k))**2 > epsmch) CYCLE
rdiag(k) = enorm(m-j, a(jp1:,k)) 1311 1285 rdiag(k) = enorm(m-j, a(jp1:,k))
wa(k) = rdiag(k) 1312 1286 wa(k) = rdiag(k)
END DO 1313 1287 END DO
rdiag(j) = -ajnorm 1314 1288 rdiag(j) = -ajnorm
END DO 1315 1289 END DO
RETURN 1316 1290 RETURN
1317 1291
! last card of subroutine qrfac. 1318 1292 ! last card of subroutine qrfac.
1319 1293
END SUBROUTINE qrfac 1320 1294 END SUBROUTINE qrfac
1321 1295
1322 1296
1323 1297 ! qrsolv
1298 ! MINPACK Subroutine (http://www.netlib.org/minpack)
1299 ! Converted to fortran90 by Alan Miller amiller @ bigpond.net.au
SUBROUTINE qrsolv(n, r, ipvt, diag, qtb, x, sdiag) 1324 1300 SUBROUTINE qrsolv(n, r, ipvt, diag, qtb, x, sdiag)
1325 1301
! N.B. Arguments LDR & WA have been removed. 1326 1302 ! N.B. Arguments LDR & WA have been removed.
INTEGER, PARAMETER :: dp = 8 1327 1303 INTEGER, PARAMETER :: dp = dp_kind
INTEGER, INTENT(IN) :: n 1328 1304 INTEGER, INTENT(IN) :: n
REAL (dp), INTENT(IN OUT) :: r(:,:) 1329 1305 REAL (dp), INTENT(IN OUT) :: r(:,:)
INTEGER, INTENT(IN) :: ipvt(:) 1330 1306 INTEGER, INTENT(IN) :: ipvt(:)
REAL (dp), INTENT(IN) :: diag(:) 1331 1307 REAL (dp), INTENT(IN) :: diag(:)
REAL (dp), INTENT(IN) :: qtb(:) 1332 1308 REAL (dp), INTENT(IN) :: qtb(:)
REAL (dp), INTENT(OUT) :: x(:) 1333 1309 REAL (dp), INTENT(OUT) :: x(:)
REAL (dp), INTENT(OUT) :: sdiag(:) 1334 1310 REAL (dp), INTENT(OUT) :: sdiag(:)
1335 1311
! ********** 1336 1312 ! **********
1337 1313
! subroutine qrsolv 1338 1314 ! subroutine qrsolv
1339 1315
! Given an m by n matrix a, an n by n diagonal matrix d, and an m-vector b, 1340 1316 ! Given an m by n matrix a, an n by n diagonal matrix d, and an m-vector b,
! the problem is to determine an x which solves the system 1341 1317 ! the problem is to determine an x which solves the system
1342 1318
! a*x = b , d*x = 0 , 1343 1319 ! a*x = b , d*x = 0 ,
1344 1320
! in the least squares sense. 1345 1321 ! in the least squares sense.
1346 1322
! This subroutine completes the solution of the problem if it is provided 1347 1323 ! This subroutine completes the solution of the problem if it is provided
! with the necessary information from the qr factorization, with column 1348 1324 ! with the necessary information from the qr factorization, with column
! pivoting, of a. That is, if a*p = q*r, where p is a permutation matrix, 1349 1325 ! pivoting, of a. That is, if a*p = q*r, where p is a permutation matrix,
! q has orthogonal columns, and r is an upper triangular matrix with diagonal 1350 1326 ! q has orthogonal columns, and r is an upper triangular matrix with diagonal
! elements of nonincreasing magnitude, then qrsolv expects the full upper 1351 1327 ! elements of nonincreasing magnitude, then qrsolv expects the full upper
! triangle of r, the permutation matrix p, and the first n components of 1352 1328 ! triangle of r, the permutation matrix p, and the first n components of
! (q transpose)*b. The system a*x = b, d*x = 0, is then equivalent to 1353 1329 ! (q transpose)*b. The system a*x = b, d*x = 0, is then equivalent to
1354 1330
! t t 1355 1331 ! t t
! r*z = q *b , p *d*p*z = 0 , 1356 1332 ! r*z = q *b , p *d*p*z = 0 ,
1357 1333
! where x = p*z. if this system does not have full rank, 1358 1334 ! where x = p*z. if this system does not have full rank,