SUEWS API Site
Documentation of SUEWS source code
Functions/Subroutines
suews_util_minpack.f95 File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine chkder (m, n, x, fvec, fjac, ldfjac, xp, fvecp, mode, err)
 
subroutine dogleg (n, r, lr, diag, qtb, delta, x)
 
real(kind=8) function enorm (n, x)
 
real(kind=8) function enorm2 (n, x)
 
subroutine fdjac1 (fcn, n, x, fvec, fjac, ldfjac, iflag, ml, mu, epsfcn, m, prms)
 
subroutine fdjac2 (fcn, m, n, x, xdat, ydat, fvec, fjac, ldfjac, iflag, epsfcn)
 
subroutine hybrd (fcn, n, x, fvec, xtol, maxfev, ml, mu, epsfcn, diag, mode, factor, nprint, info, nfev, fjac, ldfjac, r, lr, qtf, m, prms)
 
subroutine hybrd1 (fcn, n, x, fvec, tol, info, m, prms)
 
subroutine hybrj (fcn, n, x, fvec, fjac, ldfjac, xtol, maxfev, diag, mode, factor, nprint, info, nfev, njev, r, lr, qtf)
 
subroutine hybrj1 (fcn, n, x, fvec, fjac, ldfjac, tol, info)
 
subroutine lmder (fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf)
 
subroutine lmder1 (fcn, m, n, x, fvec, fjac, ldfjac, tol, info)
 
subroutine lmdif (fcn, m, n, x, xdat, ydat, fvec, ftol, xtol, gtol, maxfev, epsfcn, diag, mode, factor, nprint, info, nfev, fjac, ldfjac, ipvt, qtf)
 
subroutine lmdif1 (fcn, m, n, x, xdat, ydat, fvec, tol, info)
 
subroutine lmpar (n, r, ldr, ipvt, diag, qtb, delta, par, x, sdiag)
 
subroutine lmstr (fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf)
 
subroutine lmstr1 (fcn, m, n, x, fvec, fjac, ldfjac, tol, info)
 
subroutine qform (m, n, q, ldq)
 
subroutine qrfac (m, n, a, lda, pivot, ipvt, lipvt, rdiag, acnorm)
 
subroutine qrsolv (n, r, ldr, ipvt, diag, qtb, x, sdiag)
 
subroutine r1mpyq (m, n, a, lda, v, w)
 
subroutine r1updt (m, n, s, ls, u, v, w, sing)
 
subroutine r8vec_print (n, a, title)
 
subroutine rwupdt (n, r, ldr, w, b, alpha, c, s)
 
subroutine timestamp ()
 

Function/Subroutine Documentation

◆ chkder()

subroutine chkder ( integer(kind=4)  m,
integer(kind=4)  n,
real(kind=8), dimension(n)  x,
real(kind=8), dimension(m)  fvec,
real(kind=8), dimension(ldfjac, n)  fjac,
integer(kind=4)  ldfjac,
real(kind=8), dimension(n)  xp,
real(kind=8), dimension(m)  fvecp,
integer(kind=4)  mode,
real(kind=8), dimension(m)  err 
)

Definition at line 2 of file suews_util_minpack.f95.

2 
3 !*****************************************************************************80
4 !
5 !! CHKDER checks the gradients of M functions of N variables.
6 !
7 ! Discussion:
8 !
9 ! CHKDER checks the gradients of M nonlinear functions in N variables,
10 ! evaluated at a point X, for consistency with the functions themselves.
11 !
12 ! The user calls CHKDER twice, first with MODE = 1 and then with MODE = 2.
13 !
14 ! MODE = 1.
15 ! On input,
16 ! X contains the point of evaluation.
17 ! On output,
18 ! XP is set to a neighboring point.
19 !
20 ! Now the user must evaluate the function and gradients at X, and the
21 ! function at XP. Then the subroutine is called again:
22 !
23 ! MODE = 2.
24 ! On input,
25 ! FVEC contains the function values at X,
26 ! FJAC contains the function gradients at X.
27 ! FVECP contains the functions evaluated at XP.
28 ! On output,
29 ! ERR contains measures of correctness of the respective gradients.
30 !
31 ! The subroutine does not perform reliably if cancellation or
32 ! rounding errors cause a severe loss of significance in the
33 ! evaluation of a function. Therefore, none of the components
34 ! of X should be unusually small (in particular, zero) or any
35 ! other value which may cause loss of significance.
36 !
37 ! Licensing:
38 !
39 ! This code is distributed under the GNU LGPL license.
40 !
41 ! Modified:
42 !
43 ! 06 April 2010
44 !
45 ! Author:
46 !
47 ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom.
48 ! FORTRAN90 version by John Burkardt.
49 !
50 ! Reference:
51 !
52 ! Jorge More, Burton Garbow, Kenneth Hillstrom,
53 ! User Guide for MINPACK-1,
54 ! Technical Report ANL-80-74,
55 ! Argonne National Laboratory, 1980.
56 !
57 ! Parameters:
58 !
59 ! Input, integer ( kind = 4 ) M, is the number of functions.
60 !
61 ! Input, integer ( kind = 4 ) N, is the number of variables.
62 !
63 ! Input, real ( kind = 8 ) X(N), the point at which the jacobian is to be
64 ! evaluated.
65 !
66 ! Input, real ( kind = 8 ) FVEC(M), is used only when MODE = 2.
67 ! In that case, it should contain the function values at X.
68 !
69 ! Input, real ( kind = 8 ) FJAC(LDFJAC,N), an M by N array. When MODE = 2,
70 ! FJAC(I,J) should contain the value of dF(I)/dX(J).
71 !
72 ! Input, integer ( kind = 4 ) LDFJAC, the leading dimension of FJAC.
73 ! LDFJAC must be at least M.
74 !
75 ! Output, real ( kind = 8 ) XP(N), on output with MODE = 1, is a neighboring
76 ! point of X, at which the function is to be evaluated.
77 !
78 ! Input, real ( kind = 8 ) FVECP(M), on input with MODE = 2, is the function
79 ! value at XP.
80 !
81 ! Input, integer ( kind = 4 ) MODE, should be set to 1 on the first call, and
82 ! 2 on the second.
83 !
84 ! Output, real ( kind = 8 ) ERR(M). On output when MODE = 2, ERR contains
85 ! measures of correctness of the respective gradients. If there is no
86 ! severe loss of significance, then if ERR(I):
87 ! = 1.0D+00, the I-th gradient is correct,
88 ! = 0.0D+00, the I-th gradient is incorrect.
89 ! > 0.5D+00, the I-th gradient is probably correct.
90 ! < 0.5D+00, the I-th gradient is probably incorrect.
91 !
92  implicit none
93 
94  integer(kind=4) ldfjac
95  integer(kind=4) m
96  integer(kind=4) n
97 
98  real(kind=8) eps
99  real(kind=8) epsf
100  real(kind=8) epslog
101  real(kind=8) epsmch
102  real(kind=8) err(m)
103  real(kind=8) fjac(ldfjac, n)
104  real(kind=8) fvec(m)
105  real(kind=8) fvecp(m)
106  integer(kind=4) i
107  integer(kind=4) j
108  integer(kind=4) mode
109  real(kind=8) temp
110  real(kind=8) x(n)
111  real(kind=8) xp(n)
112 
113  epsmch = epsilon(epsmch)
114  eps = sqrt(epsmch)
115 !
116 ! MODE = 1.
117 !
118  if (mode == 1) then
119 
120  do j = 1, n
121  temp = eps*abs(x(j))
122  if (temp == 0.0d+00) then
123  temp = eps
124  end if
125  xp(j) = x(j) + temp
126  end do
127 !
128 ! MODE = 2.
129 !
130  else if (mode == 2) then
131 
132  epsf = 100.0d+00*epsmch
133  epslog = log10(eps)
134 
135  err = 0.0d+00
136 
137  do j = 1, n
138  temp = abs(x(j))
139  if (temp == 0.0d+00) then
140  temp = 1.0d+00
141  end if
142  err(1:m) = err(1:m) + temp*fjac(1:m, j)
143  end do
144 
145  do i = 1, m
146 
147  temp = 1.0d+00
148 
149  if (fvec(i) /= 0.0d+00 .and. fvecp(i) /= 0.0d+00 .and. &
150  abs(fvecp(i) - fvec(i)) >= epsf*abs(fvec(i))) then
151  temp = eps*abs((fvecp(i) - fvec(i))/eps - err(i)) &
152  /(abs(fvec(i)) + abs(fvecp(i)))
153  end if
154 
155  err(i) = 1.0d+00
156 
157  if (epsmch < temp .and. temp < eps) then
158  err(i) = (log10(temp) - epslog)/epslog
159  end if
160 
161  if (eps <= temp) then
162  err(i) = 0.0d+00
163  end if
164 
165  end do
166 
167  end if
168 
169  return

◆ dogleg()

subroutine dogleg ( integer(kind=4)  n,
real(kind=8), dimension(lr)  r,
integer(kind=4)  lr,
real(kind=8), dimension(n)  diag,
real(kind=8), dimension(n)  qtb,
real(kind=8)  delta,
real(kind=8), dimension(n)  x 
)

Definition at line 172 of file suews_util_minpack.f95.

Referenced by hybrd(), and hybrj().

172 
173 !*****************************************************************************80
174 !
175 !! DOGLEG finds the minimizing combination of Gauss-Newton and gradient steps.
176 !
177 ! Discussion:
178 !
179 ! Given an M by N matrix A, an N by N nonsingular diagonal
180 ! matrix D, an M-vector B, and a positive number DELTA, the
181 ! problem is to determine the convex combination X of the
182 ! Gauss-Newton and scaled gradient directions that minimizes
183 ! (A*X - B) in the least squares sense, subject to the
184 ! restriction that the euclidean norm of D*X be at most DELTA.
185 !
186 ! This subroutine completes the solution of the problem
187 ! if it is provided with the necessary information from the
188 ! QR factorization of A. That is, if A = Q*R, where Q has
189 ! orthogonal columns and R is an upper triangular matrix,
190 ! then DOGLEG expects the full upper triangle of R and
191 ! the first N components of Q'*B.
192 !
193 ! Licensing:
194 !
195 ! This code is distributed under the GNU LGPL license.
196 !
197 ! Modified:
198 !
199 ! 06 April 2010
200 !
201 ! Author:
202 !
203 ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom.
204 ! FORTRAN90 version by John Burkardt.
205 !
206 ! Reference:
207 !
208 ! Jorge More, Burton Garbow, Kenneth Hillstrom,
209 ! User Guide for MINPACK-1,
210 ! Technical Report ANL-80-74,
211 ! Argonne National Laboratory, 1980.
212 !
213 ! Parameters:
214 !
215 ! Input, integer ( kind = 4 ) N, the order of the matrix R.
216 !
217 ! Input, real ( kind = 8 ) R(LR), the upper triangular matrix R stored
218 ! by rows.
219 !
220 ! Input, integer ( kind = 4 ) LR, the size of the R array, which must be
221 ! no less than (N*(N+1))/2.
222 !
223 ! Input, real ( kind = 8 ) DIAG(N), the diagonal elements of the matrix D.
224 !
225 ! Input, real ( kind = 8 ) QTB(N), the first N elements of the vector Q'* B.
226 !
227 ! Input, real ( kind = 8 ) DELTA, is a positive upper bound on the
228 ! euclidean norm of D*X(1:N).
229 !
230 ! Output, real ( kind = 8 ) X(N), the desired convex combination of the
231 ! Gauss-Newton direction and the scaled gradient direction.
232 !
233  implicit none
234 
235  integer(kind=4) lr
236  integer(kind=4) n
237 
238  real(kind=8) alpha
239  real(kind=8) bnorm
240  real(kind=8) delta
241  real(kind=8) diag(n)
242  real(kind=8) enorm
243  real(kind=8) epsmch
244  real(kind=8) gnorm
245  integer(kind=4) i
246  integer(kind=4) j
247  integer(kind=4) jj
248  integer(kind=4) k
249  integer(kind=4) l
250  real(kind=8) qnorm
251  real(kind=8) qtb(n)
252  real(kind=8) r(lr)
253  real(kind=8) sgnorm
254  real(kind=8) sum2
255  real(kind=8) temp
256  real(kind=8) wa1(n)
257  real(kind=8) wa2(n)
258  real(kind=8) x(n)
259 
260  epsmch = epsilon(epsmch)
261 !
262 ! Calculate the Gauss-Newton direction.
263 !
264  jj = (n*(n + 1))/2 + 1
265 
266  do k = 1, n
267 
268  j = n - k + 1
269  jj = jj - k
270  l = jj + 1
271  sum2 = 0.0d+00
272 
273  do i = j + 1, n
274  sum2 = sum2 + r(l)*x(i)
275  l = l + 1
276  end do
277 
278  temp = r(jj)
279 
280  if (temp == 0.0d+00) then
281 
282  l = j
283  do i = 1, j
284  temp = max(temp, abs(r(l)))
285  l = l + n - i
286  end do
287 
288  if (temp == 0.0d+00) then
289  temp = epsmch
290  else
291  temp = epsmch*temp
292  end if
293 
294  end if
295 
296  x(j) = (qtb(j) - sum2)/temp
297 
298  end do
299 !
300 ! Test whether the Gauss-Newton direction is acceptable.
301 !
302  wa1(1:n) = 0.0d+00
303  wa2(1:n) = diag(1:n)*x(1:n)
304  qnorm = enorm(n, wa2)
305 
306  if (qnorm <= delta) then
307  return
308  end if
309 !
310 ! The Gauss-Newton direction is not acceptable.
311 ! Calculate the scaled gradient direction.
312 !
313  l = 1
314  do j = 1, n
315  temp = qtb(j)
316  do i = j, n
317  wa1(i) = wa1(i) + r(l)*temp
318  l = l + 1
319  end do
320  wa1(j) = wa1(j)/diag(j)
321  end do
322 !
323 ! Calculate the norm of the scaled gradient.
324 ! Test for the special case in which the scaled gradient is zero.
325 !
326  gnorm = enorm(n, wa1)
327  sgnorm = 0.0d+00
328  alpha = delta/qnorm
329 
330  if (gnorm /= 0.0d+00) then
331 !
332 ! Calculate the point along the scaled gradient which minimizes the quadratic.
333 !
334  wa1(1:n) = (wa1(1:n)/gnorm)/diag(1:n)
335 
336  l = 1
337  do j = 1, n
338  sum2 = 0.0d+00
339  do i = j, n
340  sum2 = sum2 + r(l)*wa1(i)
341  l = l + 1
342  end do
343  wa2(j) = sum2
344  end do
345 
346  temp = enorm(n, wa2)
347  sgnorm = (gnorm/temp)/temp
348 !
349 ! Test whether the scaled gradient direction is acceptable.
350 !
351  alpha = 0.0d+00
352 !
353 ! The scaled gradient direction is not acceptable.
354 ! Calculate the point along the dogleg at which the quadratic is minimized.
355 !
356  if (sgnorm < delta) then
357 
358  bnorm = enorm(n, qtb)
359  temp = (bnorm/gnorm)*(bnorm/qnorm)*(sgnorm/delta)
360  temp = temp - (delta/qnorm)*(sgnorm/delta)**2 &
361  + sqrt((temp - (delta/qnorm))**2 &
362  + (1.0d+00 - (delta/qnorm)**2) &
363  *(1.0d+00 - (sgnorm/delta)**2))
364 
365  alpha = ((delta/qnorm)*(1.0d+00 - (sgnorm/delta)**2)) &
366  /temp
367 
368  end if
369 
370  end if
371 !
372 ! Form appropriate convex combination of the Gauss-Newton
373 ! direction and the scaled gradient direction.
374 !
375  temp = (1.0d+00 - alpha)*min(sgnorm, delta)
376 
377  x(1:n) = temp*wa1(1:n) + alpha*x(1:n)
378 
379  return
real(kind=8) function enorm(n, x)
Here is the caller graph for this function:

◆ enorm()

real(kind=8) function enorm ( integer(kind=4)  n,
real(kind=8), dimension(n)  x 
)

Definition at line 382 of file suews_util_minpack.f95.

382 
383 !*****************************************************************************80
384 !
385 !! ENORM computes the Euclidean norm of a vector.
386 !
387 ! Discussion:
388 !
389 ! This is an extremely simplified version of the original ENORM
390 ! routine, which has been renamed to "ENORM2".
391 !
392 ! Licensing:
393 !
394 ! This code is distributed under the GNU LGPL license.
395 !
396 ! Modified:
397 !
398 ! 06 April 2010
399 !
400 ! Author:
401 !
402 ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom.
403 ! FORTRAN90 version by John Burkardt.
404 !
405 ! Reference:
406 !
407 ! Jorge More, Burton Garbow, Kenneth Hillstrom,
408 ! User Guide for MINPACK-1,
409 ! Technical Report ANL-80-74,
410 ! Argonne National Laboratory, 1980.
411 !
412 ! Parameters:
413 !
414 ! Input, integer ( kind = 4 ) N, is the length of the vector.
415 !
416 ! Input, real ( kind = 8 ) X(N), the vector whose norm is desired.
417 !
418 ! Output, real ( kind = 8 ) ENORM, the Euclidean norm of the vector.
419 !
420  implicit none
421 
422  integer(kind=4) n
423  real(kind=8) x(n)
424  real(kind=8) enorm
425 
426  enorm = sqrt(sum(x(1:n)**2))
427 
428  return
real(kind=8) function enorm(n, x)

◆ enorm2()

real(kind=8) function enorm2 ( integer(kind=4)  n,
real(kind=8), dimension(n)  x 
)

Definition at line 431 of file suews_util_minpack.f95.

431 
432 !*****************************************************************************80
433 !
434 !! ENORM2 computes the Euclidean norm of a vector.
435 !
436 ! Discussion:
437 !
438 ! This routine was named ENORM. It has been renamed "ENORM2",
439 ! and a simplified routine has been substituted.
440 !
441 ! The Euclidean norm is computed by accumulating the sum of
442 ! squares in three different sums. The sums of squares for the
443 ! small and large components are scaled so that no overflows
444 ! occur. Non-destructive underflows are permitted. Underflows
445 ! and overflows do not occur in the computation of the unscaled
446 ! sum of squares for the intermediate components.
447 !
448 ! The definitions of small, intermediate and large components
449 ! depend on two constants, RDWARF and RGIANT. The main
450 ! restrictions on these constants are that RDWARF^2 not
451 ! underflow and RGIANT^2 not overflow.
452 !
453 ! Licensing:
454 !
455 ! This code is distributed under the GNU LGPL license.
456 !
457 ! Modified:
458 !
459 ! 06 April 2010
460 !
461 ! Author:
462 !
463 ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom.
464 ! FORTRAN90 version by John Burkardt.
465 !
466 ! Reference:
467 !
468 ! Jorge More, Burton Garbow, Kenneth Hillstrom,
469 ! User Guide for MINPACK-1
470 ! Argonne National Laboratory,
471 ! Argonne, Illinois.
472 !
473 ! Parameters:
474 !
475 ! Input, integer ( kind = 4 ) N, is the length of the vector.
476 !
477 ! Input, real ( kind = 8 ) X(N), the vector whose norm is desired.
478 !
479 ! Output, real ( kind = 8 ) ENORM2, the Euclidean norm of the vector.
480 !
481  implicit none
482 
483  integer(kind=4) n
484 
485  real(kind=8) agiant
486  real(kind=8) enorm2
487  integer(kind=4) i
488  real(kind=8) rdwarf
489  real(kind=8) rgiant
490  real(kind=8) s1
491  real(kind=8) s2
492  real(kind=8) s3
493  real(kind=8) x(n)
494  real(kind=8) xabs
495  real(kind=8) x1max
496  real(kind=8) x3max
497 
498  rdwarf = sqrt(tiny(rdwarf))
499  rgiant = sqrt(huge(rgiant))
500 
501  s1 = 0.0d+00
502  s2 = 0.0d+00
503  s3 = 0.0d+00
504  x1max = 0.0d+00
505  x3max = 0.0d+00
506  agiant = rgiant/real(n, kind=8)
507 
508  do i = 1, n
509 
510  xabs = abs(x(i))
511 
512  if (xabs <= rdwarf) then
513 
514  if (x3max < xabs) then
515  s3 = 1.0d+00 + s3*(x3max/xabs)**2
516  x3max = xabs
517  else if (xabs /= 0.0d+00) then
518  s3 = s3 + (xabs/x3max)**2
519  end if
520 
521  else if (agiant <= xabs) then
522 
523  if (x1max < xabs) then
524  s1 = 1.0d+00 + s1*(x1max/xabs)**2
525  x1max = xabs
526  else
527  s1 = s1 + (xabs/x1max)**2
528  end if
529 
530  else
531 
532  s2 = s2 + xabs**2
533 
534  end if
535 
536  end do
537 !
538 ! Calculation of norm.
539 !
540  if (s1 /= 0.0d+00) then
541 
542  enorm2 = x1max*sqrt(s1 + (s2/x1max)/x1max)
543 
544  else if (s2 /= 0.0d+00) then
545 
546  if (x3max <= s2) then
547  enorm2 = sqrt(s2*(1.0d+00 + (x3max/s2)*(x3max*s3)))
548  else
549  enorm2 = sqrt(x3max*((s2/x3max) + (x3max*s3)))
550  end if
551 
552  else
553 
554  enorm2 = x3max*sqrt(s3)
555 
556  end if
557 
558  return
real(kind(1d0)) s2
real(kind=8) function enorm2(n, x)
real(kind(1d0)) s1

◆ fdjac1()

subroutine fdjac1 ( external  fcn,
integer(kind=4)  n,
real(kind=8), dimension(n)  x,
real(kind=8), dimension(n)  fvec,
real(kind=8), dimension(ldfjac, n)  fjac,
integer(kind=4)  ldfjac,
integer(kind=4)  iflag,
integer(kind=4)  ml,
integer(kind=4)  mu,
real(kind=8)  epsfcn,
integer(kind=4)  m,
real(kind=8), dimension(m)  prms 
)

Definition at line 561 of file suews_util_minpack.f95.

Referenced by hybrd().

561 
562 !*****************************************************************************80
563 !
564 !! FDJAC1 estimates an N by N jacobian matrix using forward differences.
565 !
566 ! Discussion:
567 !
568 ! This subroutine computes a forward-difference approximation
569 ! to the N by N jacobian matrix associated with a specified
570 ! problem of N functions in N variables. If the jacobian has
571 ! a banded form, then function evaluations are saved by only
572 ! approximating the nonzero terms.
573 !
574 ! Licensing:
575 !
576 ! This code is distributed under the GNU LGPL license.
577 !
578 ! Modified:
579 !
580 ! 06 April 2010
581 ! 30 Aug 2017: added `prms` to pass Parameters for constructing `FCN`
582 !
583 ! Author:
584 !
585 ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom.
586 ! FORTRAN90 version by John Burkardt.
587 !
588 ! Reference:
589 !
590 ! Jorge More, Burton Garbow, Kenneth Hillstrom,
591 ! User Guide for MINPACK-1,
592 ! Technical Report ANL-80-74,
593 ! Argonne National Laboratory, 1980.
594 !
595 ! Parameters:
596 !
597 ! Input, external FCN, the name of the user-supplied subroutine which
598 ! calculates the functions. The routine should have the form:
599 ! subroutine fcn ( n, x, fvec, iflag )
600 ! integer ( kind = 4 ) n
601 ! real ( kind = 8 ) fvec(n)
602 ! integer ( kind = 4 ) iflag
603 ! real ( kind = 8 ) x(n)
604 ! If IFLAG = 0 on input, then FCN is only being called to allow the user
605 ! to print out the current iterate.
606 ! If IFLAG = 1 on input, FCN should calculate the functions at X and
607 ! return this vector in FVEC.
608 ! To terminate the algorithm, FCN may set IFLAG negative on return.
609 !
610 ! Input, integer ( kind = 4 ) N, the number of functions and variables.
611 !
612 ! Input, real ( kind = 8 ) X(N), the point where the jacobian is evaluated.
613 !
614 ! Input, real ( kind = 8 ) FVEC(N), the functions evaluated at X.
615 !
616 ! Output, real ( kind = 8 ) FJAC(LDFJAC,N), the N by N approximate
617 ! jacobian matrix.
618 !
619 ! Input, integer ( kind = 4 ) LDFJAC, the leading dimension of FJAC, which
620 ! must not be less than N.
621 !
622 ! Output, integer ( kind = 4 ) IFLAG, is an error flag returned by FCN.
623 ! If FCN returns a nonzero value of IFLAG, then this routine returns
624 ! immediately to the calling program, with the value of IFLAG.
625 !
626 ! Input, integer ( kind = 4 ) ML, MU, specify the number of subdiagonals and
627 ! superdiagonals within the band of the jacobian matrix. If the
628 ! jacobian is not banded, set ML and MU to N-1.
629 !
630 ! Input, real ( kind = 8 ) EPSFCN, is used in determining a suitable step
631 ! length for the forward-difference approximation. This approximation
632 ! assumes that the relative errors in the functions are of the order of
633 ! EPSFCN. If EPSFCN is less than the machine precision, it is assumed that
634 ! the relative errors in the functions are of the order of the machine
635 ! precision.
636 !
637  implicit none
638 
639  integer(kind=4) ldfjac
640  integer(kind=4) n
641  integer(kind=4) m
642 
643  real(kind=8) eps
644  real(kind=8) epsfcn
645  real(kind=8) epsmch
646  external fcn
647  real(kind=8) fjac(ldfjac, n)
648  real(kind=8) fvec(n)
649  real(kind=8) h
650  integer(kind=4) i
651  integer(kind=4) iflag
652  integer(kind=4) j
653  integer(kind=4) k
654  integer(kind=4) ml
655  integer(kind=4) msum
656  integer(kind=4) mu
657  real(kind=8) temp
658  real(kind=8) wa1(n)
659  real(kind=8) wa2(n)
660  real(kind=8) x(n)
661  real(kind=8) prms(m)
662 
663  epsmch = epsilon(epsmch)
664 
665  eps = sqrt(max(epsfcn, epsmch))
666  msum = ml + mu + 1
667 !
668 ! Computation of dense approximate jacobian.
669 !
670  if (n <= msum) then
671 
672  do j = 1, n
673 
674  temp = x(j)
675  h = eps*abs(temp)
676  if (h == 0.0d+00) then
677  h = eps
678  end if
679 
680  iflag = 1
681  x(j) = temp + h
682  call fcn(n, x, wa1, iflag, m, prms)
683 
684  if (iflag < 0) then
685  exit
686  end if
687 
688  x(j) = temp
689  fjac(1:n, j) = (wa1(1:n) - fvec(1:n))/h
690 
691  end do
692 
693  else
694 !
695 ! Computation of banded approximate jacobian.
696 !
697  do k = 1, msum
698 
699  do j = k, n, msum
700  wa2(j) = x(j)
701  h = eps*abs(wa2(j))
702  if (h == 0.0d+00) then
703  h = eps
704  end if
705  x(j) = wa2(j) + h
706  end do
707 
708  iflag = 1
709  call fcn(n, x, wa1, iflag, m, prms)
710 
711  if (iflag < 0) then
712  exit
713  end if
714 
715  do j = k, n, msum
716 
717  x(j) = wa2(j)
718 
719  h = eps*abs(wa2(j))
720  if (h == 0.0d+00) then
721  h = eps
722  end if
723 
724  fjac(1:n, j) = 0.0d+00
725 
726  do i = 1, n
727  if (j - mu <= i .and. i <= j + ml) then
728  fjac(i, j) = (wa1(i) - fvec(i))/h
729  end if
730  end do
731 
732  end do
733 
734  end do
735 
736  end if
737 
738  return
real(kind(1d0)) h
Here is the caller graph for this function:

◆ fdjac2()

subroutine fdjac2 ( external  fcn,
integer(kind=4)  m,
integer(kind=4)  n,
real(kind=8), dimension(n)  x,
real(kind=8), dimension(m)  xdat,
real(kind=8), dimension(m)  ydat,
real(kind=8), dimension(m)  fvec,
real(kind=8), dimension(ldfjac, n)  fjac,
integer(kind=4)  ldfjac,
integer(kind=4)  iflag,
real(kind=8)  epsfcn 
)

Definition at line 741 of file suews_util_minpack.f95.

Referenced by lmdif().

741 
742 !*****************************************************************************80
743 !
744 !! FDJAC2 estimates an M by N jacobian matrix using forward differences.
745 !
746 ! Discussion:
747 !
748 ! This subroutine computes a forward-difference approximation
749 ! to the M by N jacobian matrix associated with a specified
750 ! problem of M functions in N variables.
751 !
752 ! Licensing:
753 !
754 ! This code is distributed under the GNU LGPL license.
755 !
756 ! Modified:
757 !
758 ! 06 April 2010
759 !
760 ! Author:
761 !
762 ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom.
763 ! FORTRAN90 version by John Burkardt.
764 !
765 ! Reference:
766 !
767 ! Jorge More, Burton Garbow, Kenneth Hillstrom,
768 ! User Guide for MINPACK-1,
769 ! Technical Report ANL-80-74,
770 ! Argonne National Laboratory, 1980.
771 !
772 ! Parameters:
773 !
774 ! Input, external FCN, the name of the user-supplied subroutine which
775 ! calculates the functions. The routine should have the form:
776 ! subroutine fcn ( m, n, x, xdat, ydat, fvec, iflag ) ! xdat/ydat added for AnOHM, TS 20 Jul 2017
777 ! integer ( kind = 4 ) n
778 ! real ( kind = 8 ) fvec(m)
779 ! integer ( kind = 4 ) iflag
780 ! real ( kind = 8 ) x(n)
781 !
782 ! If IFLAG = 0 on input, then FCN is only being called to allow the user
783 ! to print out the current iterate.
784 ! If IFLAG = 1 on input, FCN should calculate the functions at X and
785 ! return this vector in FVEC.
786 ! To terminate the algorithm, FCN may set IFLAG negative on return.
787 !
788 ! Input, integer ( kind = 4 ) M, is the number of functions.
789 !
790 ! Input, integer ( kind = 4 ) N, is the number of variables.
791 ! N must not exceed M.
792 !
793 ! Input, real ( kind = 8 ) X(N), the point where the jacobian is evaluated.
794 !
795 ! Input, real ( kind = 8 ) FVEC(M), the functions evaluated at X.
796 !
797 ! Output, real ( kind = 8 ) FJAC(LDFJAC,N), the M by N approximate
798 ! jacobian matrix.
799 !
800 ! Input, integer ( kind = 4 ) LDFJAC, the leading dimension of FJAC,
801 ! which must not be less than M.
802 !
803 ! Output, integer ( kind = 4 ) IFLAG, is an error flag returned by FCN.
804 ! If FCN returns a nonzero value of IFLAG, then this routine returns
805 ! immediately to the calling program, with the value of IFLAG.
806 !
807 ! Input, real ( kind = 8 ) EPSFCN, is used in determining a suitable
808 ! step length for the forward-difference approximation. This approximation
809 ! assumes that the relative errors in the functions are of the order of
810 ! EPSFCN. If EPSFCN is less than the machine precision, it is assumed that
811 ! the relative errors in the functions are of the order of the machine
812 ! precision.
813 !
814  implicit none
815 
816  integer(kind=4) ldfjac
817  integer(kind=4) m
818  integer(kind=4) n
819 
820  real(kind=8) eps
821  real(kind=8) epsfcn
822  real(kind=8) epsmch
823  external fcn
824  real(kind=8) fjac(ldfjac, n)
825  real(kind=8) fvec(m), xdat(m), ydat(m)
826  real(kind=8) h
827  ! integer ( kind = 4 ) i
828  integer(kind=4) iflag
829  integer(kind=4) j
830  real(kind=8) temp
831  real(kind=8) wa(m)
832  real(kind=8) x(n)
833 
834  epsmch = epsilon(epsmch)
835 
836  eps = sqrt(max(epsfcn, epsmch))
837 
838  do j = 1, n
839 
840  temp = x(j)
841  h = eps*abs(temp)
842  if (h == 0.0d+00) then
843  h = eps
844  end if
845 
846  iflag = 1
847  x(j) = temp + h
848  call fcn(m, n, x, xdat, ydat, wa, iflag)
849 
850  if (iflag < 0) then
851  exit
852  end if
853 
854  x(j) = temp
855  fjac(1:m, j) = (wa(1:m) - fvec(1:m))/h
856 
857  end do
858 
859  return
real(kind(1d0)) h
Here is the caller graph for this function:

◆ hybrd()

subroutine hybrd ( external  fcn,
integer(kind=4)  n,
real(kind=8), dimension(n)  x,
real(kind=8), dimension(n)  fvec,
real(kind=8)  xtol,
integer(kind=4)  maxfev,
integer(kind=4)  ml,
integer(kind=4)  mu,
real(kind=8)  epsfcn,
real(kind=8), dimension(n)  diag,
integer(kind=4)  mode,
real(kind=8)  factor,
integer(kind=4)  nprint,
integer(kind=4)  info,
integer(kind=4)  nfev,
real(kind=8), dimension(ldfjac, n)  fjac,
integer(kind=4)  ldfjac,
real(kind=8), dimension(lr)  r,
integer(kind=4)  lr,
real(kind=8), dimension(n)  qtf,
integer(kind=4)  m,
real(kind=8), dimension(m)  prms 
)

Definition at line 864 of file suews_util_minpack.f95.

References dogleg(), fdjac1(), qform(), qrfac(), r1mpyq(), and r1updt().

Referenced by hybrd1().

864 
865 !*****************************************************************************80
866 !
867 !! HYBRD seeks a zero of N nonlinear equations in N variables.
868 !
869 ! Discussion:
870 !
871 ! HYBRD finds a zero of a system of N nonlinear functions in N variables
872 ! by a modification of the Powell hybrid method. The user must provide a
873 ! subroutine which calculates the functions. The jacobian is
874 ! then calculated by a forward-difference approximation.
875 !
876 ! Licensing:
877 !
878 ! This code is distributed under the GNU LGPL license.
879 !
880 ! Modified:
881 !
882 ! 06 April 2010
883 !
884 ! Author:
885 !
886 ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom.
887 ! FORTRAN90 version by John Burkardt.
888 !
889 ! Reference:
890 !
891 ! Jorge More, Burton Garbow, Kenneth Hillstrom,
892 ! User Guide for MINPACK-1,
893 ! Technical Report ANL-80-74,
894 ! Argonne National Laboratory, 1980.
895 !
896 ! Parameters:
897 !
898 ! Input, external FCN, the name of the user-supplied subroutine which
899 ! calculates the functions. The routine should have the form:
900 ! subroutine fcn ( n, x, fvec, iflag )
901 ! integer ( kind = 4 ) n
902 ! real ( kind = 8 ) fvec(n)
903 ! integer ( kind = 4 ) iflag
904 ! real ( kind = 8 ) x(n)
905 !
906 ! If IFLAG = 0 on input, then FCN is only being called to allow the user
907 ! to print out the current iterate.
908 ! If IFLAG = 1 on input, FCN should calculate the functions at X and
909 ! return this vector in FVEC.
910 ! To terminate the algorithm, FCN may set IFLAG negative on return.
911 !
912 ! Input, integer ( kind = 4 ) N, the number of functions and variables.
913 !
914 ! Input/output, real ( kind = 8 ) X(N). On input, X must contain an initial
915 ! estimate of the solution vector. On output X contains the final
916 ! estimate of the solution vector.
917 !
918 ! Output, real ( kind = 8 ) FVEC(N), the functions evaluated at the output X.
919 !
920 ! Input, real ( kind = 8 ) XTOL. Termination occurs when the relative error
921 ! between two consecutive iterates is at most XTOL. XTOL should be
922 ! nonnegative.
923 !
924 ! Input, integer ( kind = 4 ) MAXFEV. Termination occurs when the number of
925 ! calls to FCN is at least MAXFEV by the end of an iteration.
926 !
927 ! Input, integer ( kind = 4 ) ML, MU, specify the number of subdiagonals and
928 ! superdiagonals within the band of the jacobian matrix. If the jacobian
929 ! is not banded, set ML and MU to at least n - 1.
930 !
931 ! Input, real ( kind = 8 ) EPSFCN, is used in determining a suitable step
932 ! length for the forward-difference approximation. This approximation
933 ! assumes that the relative errors in the functions are of the order of
934 ! EPSFCN. If EPSFCN is less than the machine precision, it is assumed that
935 ! the relative errors in the functions are of the order of the machine
936 ! precision.
937 !
938 ! Input/output, real ( kind = 8 ) DIAG(N). If MODE = 1, then DIAG is set
939 ! internally. If MODE = 2, then DIAG must contain positive entries that
940 ! serve as multiplicative scale factors for the variables.
941 !
942 ! Input, integer ( kind = 4 ) MODE, scaling option.
943 ! 1, variables will be scaled internally.
944 ! 2, scaling is specified by the input DIAG vector.
945 !
946 ! Input, real ( kind = 8 ) FACTOR, determines the initial step bound. This
947 ! bound is set to the product of FACTOR and the euclidean norm of DIAG*X if
948 ! nonzero, or else to FACTOR itself. In most cases, FACTOR should lie
949 ! in the interval (0.1, 100) with 100 the recommended value.
950 !
951 ! Input, integer ( kind = 4 ) NPRINT, enables controlled printing of
952 ! iterates if it is positive. In this case, FCN is called with IFLAG = 0 at
953 ! the beginning of the first iteration and every NPRINT iterations thereafter
954 ! and immediately prior to return, with X and FVEC available
955 ! for printing. If NPRINT is not positive, no special calls
956 ! of FCN with IFLAG = 0 are made.
957 !
958 ! Output, integer ( kind = 4 ) INFO, error flag. If the user has terminated
959 ! execution, INFO is set to the (negative) value of IFLAG.
960 ! See the description of FCN.
961 ! Otherwise, INFO is set as follows:
962 ! 0, improper input parameters.
963 ! 1, relative error between two consecutive iterates is at most XTOL.
964 ! 2, number of calls to FCN has reached or exceeded MAXFEV.
965 ! 3, XTOL is too small. No further improvement in the approximate
966 ! solution X is possible.
967 ! 4, iteration is not making good progress, as measured by the improvement
968 ! from the last five jacobian evaluations.
969 ! 5, iteration is not making good progress, as measured by the improvement
970 ! from the last ten iterations.
971 !
972 ! Output, integer ( kind = 4 ) NFEV, the number of calls to FCN.
973 !
974 ! Output, real ( kind = 8 ) FJAC(LDFJAC,N), an N by N array which contains
975 ! the orthogonal matrix Q produced by the QR factorization of the final
976 ! approximate jacobian.
977 !
978 ! Input, integer ( kind = 4 ) LDFJAC, the leading dimension of FJAC.
979 ! LDFJAC must be at least N.
980 !
981 ! Output, real ( kind = 8 ) R(LR), the upper triangular matrix produced by
982 ! the QR factorization of the final approximate jacobian, stored rowwise.
983 !
984 ! Input, integer ( kind = 4 ) LR, the size of the R array, which must be no
985 ! less than (N*(N+1))/2.
986 !
987 ! Output, real ( kind = 8 ) QTF(N), contains the vector Q'*FVEC.
988 !
989  implicit none
990 
991  integer(kind=4) ldfjac
992  integer(kind=4) lr
993  integer(kind=4) n
994  integer(kind=4) m
995 
996  real(kind=8) actred
997  real(kind=8) delta
998  real(kind=8) diag(n)
999  real(kind=8) enorm
1000  real(kind=8) epsfcn
1001  real(kind=8) epsmch
1002  real(kind=8) factor
1003  external fcn
1004  real(kind=8) fjac(ldfjac, n)
1005  real(kind=8) fnorm
1006  real(kind=8) fnorm1
1007  real(kind=8) fvec(n)
1008  integer(kind=4) i
1009  integer(kind=4) iflag
1010  integer(kind=4) info
1011  integer(kind=4) iter
1012  integer(kind=4) iwa(1)
1013  integer(kind=4) j
1014  logical jeval
1015  integer(kind=4) l
1016  integer(kind=4) maxfev
1017  integer(kind=4) ml
1018  integer(kind=4) mode
1019  integer(kind=4) msum
1020  integer(kind=4) mu
1021  integer(kind=4) ncfail
1022  integer(kind=4) nslow1
1023  integer(kind=4) nslow2
1024  integer(kind=4) ncsuc
1025  integer(kind=4) nfev
1026  integer(kind=4) nprint
1027  logical pivot
1028  real(kind=8) pnorm
1029  real(kind=8) prered
1030  real(kind=8) qtf(n)
1031  real(kind=8) r(lr)
1032  real(kind=8) ratio
1033  logical sing
1034  real(kind=8) sum2
1035  real(kind=8) temp
1036  real(kind=8) wa1(n)
1037  real(kind=8) wa2(n)
1038  real(kind=8) wa3(n)
1039  real(kind=8) wa4(n)
1040  real(kind=8) x(n)
1041  real(kind=8) xnorm
1042  real(kind=8) xtol
1043  real(kind=8) prms(m)
1044 
1045  epsmch = epsilon(epsmch)
1046 
1047  info = 0
1048  iflag = 0
1049  nfev = 0
1050 !
1051 ! Check the input parameters for errors.
1052 !
1053  if (n <= 0) then
1054  return
1055  else if (xtol < 0.0d+00) then
1056  return
1057  else if (maxfev <= 0) then
1058  return
1059  else if (ml < 0) then
1060  return
1061  else if (mu < 0) then
1062  return
1063  else if (factor <= 0.0d+00) then
1064  return
1065  else if (ldfjac < n) then
1066  return
1067  else if (lr < (n*(n + 1))/2) then
1068  return
1069  end if
1070 
1071  if (mode == 2) then
1072 
1073  do j = 1, n
1074  if (diag(j) <= 0.0d+00) then
1075  go to 300
1076  end if
1077  end do
1078 
1079  end if
1080 !
1081 ! Evaluate the function at the starting point
1082 ! and calculate its norm.
1083 !
1084  iflag = 1
1085  call fcn(n, x, fvec, iflag, m, prms)
1086  nfev = 1
1087 
1088  if (iflag < 0) then
1089  go to 300
1090  end if
1091 
1092  fnorm = enorm(n, fvec)
1093 !
1094 ! Determine the number of calls to FCN needed to compute the jacobian matrix.
1095 !
1096  msum = min(ml + mu + 1, n)
1097 !
1098 ! Initialize iteration counter and monitors.
1099 !
1100  iter = 1
1101  ncsuc = 0
1102  ncfail = 0
1103  nslow1 = 0
1104  nslow2 = 0
1105 !
1106 ! Beginning of the outer loop.
1107 !
1108 30 continue
1109 
1110  jeval = .true.
1111 !
1112 ! Calculate the jacobian matrix.
1113 !
1114  iflag = 2
1115  call fdjac1(fcn, n, x, fvec, fjac, ldfjac, iflag, ml, mu, epsfcn, m, prms)
1116 
1117  nfev = nfev + msum
1118 
1119  if (iflag < 0) then
1120  go to 300
1121  end if
1122 !
1123 ! Compute the QR factorization of the jacobian.
1124 !
1125  pivot = .false.
1126  call qrfac(n, n, fjac, ldfjac, pivot, iwa, 1, wa1, wa2)
1127 !
1128 ! On the first iteration, if MODE is 1, scale according
1129 ! to the norms of the columns of the initial jacobian.
1130 !
1131  if (iter == 1) then
1132 
1133  if (mode /= 2) then
1134 
1135  diag(1:n) = wa2(1:n)
1136  do j = 1, n
1137  if (wa2(j) == 0.0d+00) then
1138  diag(j) = 1.0d+00
1139  end if
1140  end do
1141 
1142  end if
1143 !
1144 ! On the first iteration, calculate the norm of the scaled X
1145 ! and initialize the step bound DELTA.
1146 !
1147  wa3(1:n) = diag(1:n)*x(1:n)
1148  xnorm = enorm(n, wa3)
1149  delta = factor*xnorm
1150  if (delta == 0.0d+00) then
1151  delta = factor
1152  end if
1153 
1154  end if
1155 !
1156 ! Form Q' * FVEC and store in QTF.
1157 !
1158  qtf(1:n) = fvec(1:n)
1159 
1160  do j = 1, n
1161 
1162  if (fjac(j, j) /= 0.0d+00) then
1163  temp = -dot_product(qtf(j:n), fjac(j:n, j))/fjac(j, j)
1164  qtf(j:n) = qtf(j:n) + fjac(j:n, j)*temp
1165  end if
1166 
1167  end do
1168 !
1169 ! Copy the triangular factor of the QR factorization into R.
1170 !
1171  sing = .false.
1172 
1173  do j = 1, n
1174  l = j
1175  do i = 1, j - 1
1176  r(l) = fjac(i, j)
1177  l = l + n - i
1178  end do
1179  r(l) = wa1(j)
1180  if (wa1(j) == 0.0d+00) then
1181  sing = .true.
1182  end if
1183  end do
1184 !
1185 ! Accumulate the orthogonal factor in FJAC.
1186 !
1187  call qform(n, n, fjac, ldfjac)
1188 !
1189 ! Rescale if necessary.
1190 !
1191  if (mode /= 2) then
1192  do j = 1, n
1193  diag(j) = max(diag(j), wa2(j))
1194  end do
1195  end if
1196 !
1197 ! Beginning of the inner loop.
1198 !
1199 180 continue
1200 !
1201 ! If requested, call FCN to enable printing of iterates.
1202 !
1203  if (0 < nprint) then
1204  iflag = 0
1205  if (mod(iter - 1, nprint) == 0) then
1206  call fcn(n, x, fvec, iflag, m, prms)
1207  end if
1208  if (iflag < 0) then
1209  go to 300
1210  end if
1211  end if
1212 !
1213 ! Determine the direction P.
1214 !
1215  call dogleg(n, r, lr, diag, qtf, delta, wa1)
1216 !
1217 ! Store the direction P and X + P.
1218 ! Calculate the norm of P.
1219 !
1220  wa1(1:n) = -wa1(1:n)
1221  wa2(1:n) = x(1:n) + wa1(1:n)
1222  wa3(1:n) = diag(1:n)*wa1(1:n)
1223 
1224  pnorm = enorm(n, wa3)
1225 !
1226 ! On the first iteration, adjust the initial step bound.
1227 !
1228  if (iter == 1) then
1229  delta = min(delta, pnorm)
1230  end if
1231 !
1232 ! Evaluate the function at X + P and calculate its norm.
1233 !
1234  iflag = 1
1235  call fcn(n, wa2, wa4, iflag, m, prms)
1236  nfev = nfev + 1
1237 
1238  if (iflag < 0) then
1239  go to 300
1240  end if
1241 
1242  fnorm1 = enorm(n, wa4)
1243 !
1244 ! Compute the scaled actual reduction.
1245 !
1246  actred = -1.0d+00
1247  if (fnorm1 < fnorm) then
1248  actred = 1.0d+00 - (fnorm1/fnorm)**2
1249  endif
1250 !
1251 ! Compute the scaled predicted reduction.
1252 !
1253  l = 1
1254  do i = 1, n
1255  sum2 = 0.0d+00
1256  do j = i, n
1257  sum2 = sum2 + r(l)*wa1(j)
1258  l = l + 1
1259  end do
1260  wa3(i) = qtf(i) + sum2
1261  end do
1262 
1263  temp = enorm(n, wa3)
1264  prered = 0.0d+00
1265  if (temp < fnorm) then
1266  prered = 1.0d+00 - (temp/fnorm)**2
1267  end if
1268 !
1269 ! Compute the ratio of the actual to the predicted reduction.
1270 !
1271  ratio = 0.0d+00
1272  if (0.0d+00 < prered) then
1273  ratio = actred/prered
1274  end if
1275 !
1276 ! Update the step bound.
1277 !
1278  if (ratio < 0.1d+00) then
1279 
1280  ncsuc = 0
1281  ncfail = ncfail + 1
1282  delta = 0.5d+00*delta
1283 
1284  else
1285 
1286  ncfail = 0
1287  ncsuc = ncsuc + 1
1288 
1289  if (0.5d+00 <= ratio .or. 1 < ncsuc) then
1290  delta = max(delta, pnorm/0.5d+00)
1291  end if
1292 
1293  if (abs(ratio - 1.0d+00) <= 0.1d+00) then
1294  delta = pnorm/0.5d+00
1295  end if
1296 
1297  end if
1298 !
1299 ! Test for successful iteration.
1300 !
1301 ! Successful iteration.
1302 ! Update X, FVEC, and their norms.
1303 !
1304  if (0.0001d+00 <= ratio) then
1305  x(1:n) = wa2(1:n)
1306  wa2(1:n) = diag(1:n)*x(1:n)
1307  fvec(1:n) = wa4(1:n)
1308  xnorm = enorm(n, wa2)
1309  fnorm = fnorm1
1310  iter = iter + 1
1311  end if
1312 !
1313 ! Determine the progress of the iteration.
1314 !
1315  nslow1 = nslow1 + 1
1316  if (0.001d+00 <= actred) then
1317  nslow1 = 0
1318  end if
1319 
1320  if (jeval) then
1321  nslow2 = nslow2 + 1
1322  end if
1323 
1324  if (0.1d+00 <= actred) then
1325  nslow2 = 0
1326  end if
1327 !
1328 ! Test for convergence.
1329 !
1330  if (delta <= xtol*xnorm .or. fnorm == 0.0d+00) then
1331  info = 1
1332  end if
1333 
1334  if (info /= 0) then
1335  go to 300
1336  end if
1337 !
1338 ! Tests for termination and stringent tolerances.
1339 !
1340  if (maxfev <= nfev) then
1341  info = 2
1342  end if
1343 
1344  if (0.1d+00*max(0.1d+00*delta, pnorm) <= epsmch*xnorm) then
1345  info = 3
1346  end if
1347 
1348  if (nslow2 == 5) then
1349  info = 4
1350  end if
1351 
1352  if (nslow1 == 10) then
1353  info = 5
1354  end if
1355 
1356  if (info /= 0) then
1357  go to 300
1358  end if
1359 !
1360 ! Criterion for recalculating jacobian approximation
1361 ! by forward differences.
1362 !
1363  if (ncfail == 2) then
1364  go to 290
1365  end if
1366 !
1367 ! Calculate the rank one modification to the jacobian
1368 ! and update QTF if necessary.
1369 !
1370  do j = 1, n
1371  sum2 = dot_product(wa4(1:n), fjac(1:n, j))
1372  wa2(j) = (sum2 - wa3(j))/pnorm
1373  wa1(j) = diag(j)*((diag(j)*wa1(j))/pnorm)
1374  if (0.0001d+00 <= ratio) then
1375  qtf(j) = sum2
1376  end if
1377  end do
1378 !
1379 ! Compute the QR factorization of the updated jacobian.
1380 !
1381  call r1updt(n, n, r, lr, wa1, wa2, wa3, sing)
1382  call r1mpyq(n, n, fjac, ldfjac, wa2, wa3)
1383  call r1mpyq(1, n, qtf, 1, wa2, wa3)
1384 !
1385 ! End of the inner loop.
1386 !
1387  jeval = .false.
1388  go to 180
1389 
1390 290 continue
1391 !
1392 ! End of the outer loop.
1393 !
1394  go to 30
1395 
1396 300 continue
1397 !
1398 ! Termination, either normal or user imposed.
1399 !
1400  if (iflag < 0) then
1401  info = iflag
1402  end if
1403 
1404  iflag = 0
1405 
1406  if (0 < nprint) then
1407  call fcn(n, x, fvec, iflag, m, prms)
1408  end if
1409 
1410  return
subroutine dogleg(n, r, lr, diag, qtb, delta, x)
subroutine r1updt(m, n, s, ls, u, v, w, sing)
real(kind=8) function enorm(n, x)
subroutine qrfac(m, n, a, lda, pivot, ipvt, lipvt, rdiag, acnorm)
subroutine r1mpyq(m, n, a, lda, v, w)
subroutine fdjac1(fcn, n, x, fvec, fjac, ldfjac, iflag, ml, mu, epsfcn, m, prms)
subroutine qform(m, n, q, ldq)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ hybrd1()

subroutine hybrd1 ( external  fcn,
integer(kind=4)  n,
real(kind=8), dimension(n)  x,
real(kind=8), dimension(n)  fvec,
real(kind=8)  tol,
integer(kind=4)  info,
integer(kind=4)  m,
real(kind=8), dimension(m)  prms 
)

Definition at line 1414 of file suews_util_minpack.f95.

References hybrd().

Referenced by anohm_module::anohm_bo_cal().

1414 
1415 !*****************************************************************************80
1416 !
1417 !! HYBRD1 seeks a zero of N nonlinear equations in N variables.
1418 !
1419 ! Discussion:
1420 !
1421 ! HYBRD1 finds a zero of a system of N nonlinear functions in N variables
1422 ! by a modification of the Powell hybrid method. This is done by using the
1423 ! more general nonlinear equation solver HYBRD. The user must provide a
1424 ! subroutine which calculates the functions. The jacobian is then
1425 ! calculated by a forward-difference approximation.
1426 !
1427 ! Licensing:
1428 !
1429 ! This code is distributed under the GNU LGPL license.
1430 !
1431 ! Modified:
1432 !
1433 ! 19 August 2016
1434 ! 30 Aug 2017: added `prms` to pass paramters for constructing `fcn`
1435 !
1436 ! Author:
1437 !
1438 ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom.
1439 ! FORTRAN90 version by John Burkardt.
1440 !
1441 ! Reference:
1442 !
1443 ! Jorge More, Burton Garbow, Kenneth Hillstrom,
1444 ! User Guide for MINPACK-1,
1445 ! Technical Report ANL-80-74,
1446 ! Argonne National Laboratory, 1980.
1447 !
1448 ! Parameters:
1449 !
1450 ! Input, external FCN, the name of the user-supplied subroutine which
1451 ! calculates the functions. The routine should have the form:
1452 ! subroutine fcn ( n, x, fvec, iflag )
1453 ! integer ( kind = 4 ) n
1454 ! real ( kind = 8 ) fvec(n)
1455 ! integer ( kind = 4 ) iflag
1456 ! real ( kind = 8 ) x(n)
1457 ! integer ( kind = 4 ) m
1458 ! real ( kind = 8 ) prm(m)
1459 ! If IFLAG = 0 on input, then FCN is only being called to allow the user
1460 ! to print out the current iterate.
1461 ! If IFLAG = 1 on input, FCN should calculate the functions at X and
1462 ! return this vector in FVEC.
1463 ! To terminate the algorithm, FCN may set IFLAG negative on return.
1464 !
1465 ! Input, integer ( kind = 4 ) N, the number of functions and variables.
1466 !
1467 ! Input/output, real ( kind = 8 ) X(N). On input, X must contain an initial
1468 ! estimate of the solution vector. On output X contains the final
1469 ! estimate of the solution vector.
1470 !
1471 ! Input, real ( kind = 8 ) FVEC(N), the functions evaluated at the output X.
1472 !
1473 ! Input, integer ( kind = 4 ) M, the number of parameters for constructing FCN.
1474 !
1475 ! Input, real ( kind = 8 ) PRMS(M), static paramters for constructing FCN.
1476 !
1477 ! Input, real ( kind = 8 ) TOL. Termination occurs when the algorithm
1478 ! estimates that the relative error between X and the solution is at
1479 ! most TOL. TOL should be nonnegative.
1480 !
1481 ! Output, integer ( kind = 4 ) INFO, error flag. If the user has terminated
1482 ! execution, INFO is set to the (negative) value of IFLAG. See the
1483 ! description of FCN.
1484 ! Otherwise, INFO is set as follows:
1485 ! 0, improper input parameters.
1486 ! 1, algorithm estimates that the relative error between X and the
1487 ! solution is at most TOL.
1488 ! 2, number of calls to FCN has reached or exceeded 200*(N+1).
1489 ! 3, TOL is too small. No further improvement in the approximate
1490 ! solution X is possible.
1491 ! 4, the iteration is not making good progress.
1492 !
1493  implicit none
1494 
1495  ! integer ( kind = 4 ) lwa
1496  integer(kind=4) n
1497  integer(kind=4) m
1498 
1499  real(kind=8) diag(n)
1500  real(kind=8) epsfcn
1501  real(kind=8) factor
1502  external fcn
1503  real(kind=8) fjac(n, n)
1504  real(kind=8) fvec(n)
1505  integer(kind=4) info
1506  ! integer ( kind = 4 ) j
1507  integer(kind=4) ldfjac
1508  integer(kind=4) lr
1509  integer(kind=4) maxfev
1510  integer(kind=4) ml
1511  integer(kind=4) mode
1512  integer(kind=4) mu
1513  integer(kind=4) nfev
1514  integer(kind=4) nprint
1515  real(kind=8) qtf(n)
1516  real(kind=8) r((n*(n + 1))/2)
1517  real(kind=8) tol
1518  real(kind=8) x(n)
1519  real(kind=8) xtol
1520  real(kind=8) prms(m)
1521 
1522  if (n <= 0) then
1523  info = 0
1524  return
1525  end if
1526 
1527  if (tol < 0.0d+00) then
1528  info = 0
1529  return
1530  end if
1531 
1532  xtol = tol
1533  maxfev = 200*(n + 1)
1534  ml = n - 1
1535  mu = n - 1
1536  epsfcn = 0.0d+00
1537  diag(1:n) = 1.0d+00
1538  mode = 2
1539  factor = 100.0d+00
1540  nprint = 0
1541  info = 0
1542  nfev = 0
1543  fjac(1:n, 1:n) = 0.0d+00
1544  ldfjac = n
1545  r(1:(n*(n + 1))/2) = 0.0d+00
1546  lr = (n*(n + 1))/2
1547  qtf(1:n) = 0.0d+00
1548 
1549  call hybrd(fcn, n, x, fvec, xtol, maxfev, ml, mu, epsfcn, diag, mode, &
1550  factor, nprint, info, nfev, fjac, ldfjac, r, lr, qtf, m, prms)
1551 
1552  if (info == 5) then
1553  info = 4
1554  end if
1555 
1556  return
subroutine hybrd(fcn, n, x, fvec, xtol, maxfev, ml, mu, epsfcn, diag, mode, factor, nprint, info, nfev, fjac, ldfjac, r, lr, qtf, m, prms)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ hybrj()

subroutine hybrj ( external  fcn,
integer(kind=4)  n,
real(kind=8), dimension(n)  x,
real(kind=8), dimension(n)  fvec,
real(kind=8), dimension(ldfjac, n)  fjac,
integer(kind=4)  ldfjac,
real(kind=8)  xtol,
integer(kind=4)  maxfev,
real(kind=8), dimension(n)  diag,
integer(kind=4)  mode,
real(kind=8)  factor,
integer(kind=4)  nprint,
integer(kind=4)  info,
integer(kind=4)  nfev,
integer(kind=4)  njev,
real(kind=8), dimension(lr)  r,
integer(kind=4)  lr,
real(kind=8), dimension(n)  qtf 
)

Definition at line 1560 of file suews_util_minpack.f95.

References dogleg(), qform(), qrfac(), r1mpyq(), and r1updt().

Referenced by hybrj1().

1560 
1561 !*****************************************************************************80
1562 !
1563 !! HYBRJ seeks a zero of N nonlinear equations in N variables.
1564 !
1565 ! Discussion:
1566 !
1567 ! HYBRJ finds a zero of a system of N nonlinear functions in N variables
1568 ! by a modification of the Powell hybrid method. The user must provide a
1569 ! subroutine which calculates the functions and the jacobian.
1570 !
1571 ! Licensing:
1572 !
1573 ! This code is distributed under the GNU LGPL license.
1574 !
1575 ! Modified:
1576 !
1577 ! 06 April 2010
1578 !
1579 ! Author:
1580 !
1581 ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom.
1582 ! FORTRAN90 version by John Burkardt.
1583 !
1584 ! Reference:
1585 !
1586 ! Jorge More, Burton Garbow, Kenneth Hillstrom,
1587 ! User Guide for MINPACK-1,
1588 ! Technical Report ANL-80-74,
1589 ! Argonne National Laboratory, 1980.
1590 !
1591 ! Parameters:
1592 !
1593 ! Input, external FCN, the name of the user-supplied subroutine which
1594 ! calculates the functions and the jacobian. FCN should have the form:
1595 !
1596 ! subroutine fcn ( n, x, fvec, fjac, ldfjac, iflag )
1597 ! integer ( kind = 4 ) ldfjac
1598 ! integer ( kind = 4 ) n
1599 ! real ( kind = 8 ) fjac(ldfjac,n)
1600 ! real ( kind = 8 ) fvec(n)
1601 ! integer ( kind = 4 ) iflag
1602 ! real ( kind = 8 ) x(n)
1603 !
1604 ! If IFLAG = 0 on input, then FCN is only being called to allow the user
1605 ! to print out the current iterate.
1606 ! If IFLAG = 1 on input, FCN should calculate the functions at X and
1607 ! return this vector in FVEC.
1608 ! If IFLAG = 2 on input, FCN should calculate the jacobian at X and
1609 ! return this matrix in FJAC.
1610 ! To terminate the algorithm, FCN may set IFLAG negative on return.
1611 !
1612 ! Input, integer ( kind = 4 ) N, the number of functions and variables.
1613 !
1614 ! Input/output, real ( kind = 8 ) X(N). On input, X must contain an initial
1615 ! estimate of the solution vector. On output X contains the final
1616 ! estimate of the solution vector.
1617 !
1618 ! Output, real ( kind = 8 ) FVEC(N), the functions evaluated at the output X.
1619 !
1620 ! Output, real ( kind = 8 ) FJAC(LDFJAC,N), an N by N matrix, containing
1621 ! the orthogonal matrix Q produced by the QR factorization
1622 ! of the final approximate jacobian.
1623 !
1624 ! Input, integer ( kind = 4 ) LDFJAC, the leading dimension of the
1625 ! array FJAC. LDFJAC must be at least N.
1626 !
1627 ! Input, real ( kind = 8 ) XTOL. Termination occurs when the relative error
1628 ! between two consecutive iterates is at most XTOL. XTOL should be
1629 ! nonnegative.
1630 !
1631 ! Input, integer ( kind = 4 ) MAXFEV. Termination occurs when the number of
1632 ! calls to FCN is at least MAXFEV by the end of an iteration.
1633 !
1634 ! Input/output, real ( kind = 8 ) DIAG(N). If MODE = 1, then DIAG is set
1635 ! internally. If MODE = 2, then DIAG must contain positive entries that
1636 ! serve as multiplicative scale factors for the variables.
1637 !
1638 ! Input, integer ( kind = 4 ) MODE, scaling option.
1639 ! 1, variables will be scaled internally.
1640 ! 2, scaling is specified by the input DIAG vector.
1641 !
1642 ! Input, real ( kind = 8 ) FACTOR, determines the initial step bound. This
1643 ! bound is set to the product of FACTOR and the euclidean norm of DIAG*X if
1644 ! nonzero, or else to FACTOR itself. In most cases, FACTOR should lie
1645 ! in the interval (0.1, 100) with 100 the recommended value.
1646 !
1647 ! Input, integer ( kind = 4 ) NPRINT, enables controlled printing of iterates
1648 ! if it is positive. In this case, FCN is called with IFLAG = 0 at the
1649 ! beginning of the first iteration and every NPRINT iterations thereafter
1650 ! and immediately prior to return, with X and FVEC available
1651 ! for printing. If NPRINT is not positive, no special calls
1652 ! of FCN with IFLAG = 0 are made.
1653 !
1654 ! Output, integer ( kind = 4 ) INFO, error flag. If the user has terminated
1655 ! execution, INFO is set to the (negative) value of IFLAG.
1656 ! See the description of FCN. Otherwise, INFO is set as follows:
1657 ! 0, improper input parameters.
1658 ! 1, relative error between two consecutive iterates is at most XTOL.
1659 ! 2, number of calls to FCN with IFLAG = 1 has reached MAXFEV.
1660 ! 3, XTOL is too small. No further improvement in
1661 ! the approximate solution X is possible.
1662 ! 4, iteration is not making good progress, as measured by the
1663 ! improvement from the last five jacobian evaluations.
1664 ! 5, iteration is not making good progress, as measured by the
1665 ! improvement from the last ten iterations.
1666 !
1667 ! Output, integer ( kind = 4 ) NFEV, the number of calls to FCN
1668 ! with IFLAG = 1.
1669 !
1670 ! Output, integer ( kind = 4 ) NJEV, the number of calls to FCN
1671 ! with IFLAG = 2.
1672 !
1673 ! Output, real ( kind = 8 ) R(LR), the upper triangular matrix produced
1674 ! by the QR factorization of the final approximate jacobian, stored rowwise.
1675 !
1676 ! Input, integer ( kind = 4 ) LR, the size of the R array, which must
1677 ! be no less than (N*(N+1))/2.
1678 !
1679 ! Output, real ( kind = 8 ) QTF(N), contains the vector Q'*FVEC.
1680 !
1681  implicit none
1682 
1683  integer(kind=4) ldfjac
1684  integer(kind=4) lr
1685  integer(kind=4) n
1686 
1687  real(kind=8) actred
1688  real(kind=8) delta
1689  real(kind=8) diag(n)
1690  real(kind=8) enorm
1691  real(kind=8) epsmch
1692  real(kind=8) factor
1693  external fcn
1694  real(kind=8) fjac(ldfjac, n)
1695  real(kind=8) fnorm
1696  real(kind=8) fnorm1
1697  real(kind=8) fvec(n)
1698  integer(kind=4) i
1699  integer(kind=4) iflag
1700  integer(kind=4) info
1701  integer(kind=4) iter
1702  integer(kind=4) iwa(1)
1703  integer(kind=4) j
1704  logical jeval
1705  integer(kind=4) l
1706  integer(kind=4) maxfev
1707  integer(kind=4) mode
1708  integer(kind=4) ncfail
1709  integer(kind=4) nslow1
1710  integer(kind=4) nslow2
1711  integer(kind=4) ncsuc
1712  integer(kind=4) nfev
1713  integer(kind=4) njev
1714  integer(kind=4) nprint
1715  logical pivot
1716  real(kind=8) pnorm
1717  real(kind=8) prered
1718  real(kind=8) qtf(n)
1719  real(kind=8) r(lr)
1720  real(kind=8) ratio
1721  logical sing
1722  real(kind=8) sum2
1723  real(kind=8) temp
1724  real(kind=8) wa1(n)
1725  real(kind=8) wa2(n)
1726  real(kind=8) wa3(n)
1727  real(kind=8) wa4(n)
1728  real(kind=8) x(n)
1729  real(kind=8) xnorm
1730  real(kind=8) xtol
1731 
1732  epsmch = epsilon(epsmch)
1733 
1734  info = 0
1735  iflag = 0
1736  nfev = 0
1737  njev = 0
1738 !
1739 ! Check the input parameters for errors.
1740 !
1741  if (n <= 0) then
1742  if (iflag < 0) then
1743  info = iflag
1744  end if
1745  iflag = 0
1746  if (0 < nprint) then
1747  call fcn(n, x, fvec, fjac, ldfjac, iflag)
1748  end if
1749  return
1750  end if
1751 
1752  if (ldfjac < n .or. &
1753  xtol < 0.0d+00 .or. &
1754  maxfev <= 0 .or. &
1755  factor <= 0.0d+00 .or. &
1756  lr < (n*(n + 1))/2) then
1757  if (iflag < 0) then
1758  info = iflag
1759  end if
1760  iflag = 0
1761  if (0 < nprint) then
1762  call fcn(n, x, fvec, fjac, ldfjac, iflag)
1763  end if
1764  return
1765  end if
1766 
1767  if (mode == 2) then
1768  do j = 1, n
1769  if (diag(j) <= 0.0d+00) then
1770  if (iflag < 0) then
1771  info = iflag
1772  end if
1773  iflag = 0
1774  if (0 < nprint) then
1775  call fcn(n, x, fvec, fjac, ldfjac, iflag)
1776  end if
1777  return
1778  end if
1779  end do
1780  end if
1781 !
1782 ! Evaluate the function at the starting point and calculate its norm.
1783 !
1784  iflag = 1
1785  call fcn(n, x, fvec, fjac, ldfjac, iflag)
1786  nfev = 1
1787 
1788  if (iflag < 0) then
1789  if (iflag < 0) then
1790  info = iflag
1791  end if
1792  iflag = 0
1793  if (0 < nprint) then
1794  call fcn(n, x, fvec, fjac, ldfjac, iflag)
1795  end if
1796  return
1797  end if
1798 
1799  fnorm = enorm(n, fvec)
1800 !
1801 ! Initialize iteration counter and monitors.
1802 !
1803  iter = 1
1804  ncsuc = 0
1805  ncfail = 0
1806  nslow1 = 0
1807  nslow2 = 0
1808 !
1809 ! Beginning of the outer loop.
1810 !
1811  do
1812 
1813  jeval = .true.
1814 !
1815 ! Calculate the jacobian matrix.
1816 !
1817  iflag = 2
1818  call fcn(n, x, fvec, fjac, ldfjac, iflag)
1819  njev = njev + 1
1820 
1821  if (iflag < 0) then
1822  info = iflag
1823  iflag = 0
1824  if (0 < nprint) then
1825  call fcn(n, x, fvec, fjac, ldfjac, iflag)
1826  end if
1827  return
1828  end if
1829 !
1830 ! Compute the QR factorization of the jacobian.
1831 !
1832  pivot = .false.
1833  call qrfac(n, n, fjac, ldfjac, pivot, iwa, 1, wa1, wa2)
1834 !
1835 ! On the first iteration, if MODE is 1, scale according
1836 ! to the norms of the columns of the initial jacobian.
1837 !
1838  if (iter == 1) then
1839 
1840  if (mode /= 2) then
1841  diag(1:n) = wa2(1:n)
1842  do j = 1, n
1843  if (wa2(j) == 0.0d+00) then
1844  diag(j) = 1.0d+00
1845  end if
1846  end do
1847  end if
1848 !
1849 ! On the first iteration, calculate the norm of the scaled X
1850 ! and initialize the step bound DELTA.
1851 !
1852  wa3(1:n) = diag(1:n)*x(1:n)
1853  xnorm = enorm(n, wa3)
1854  delta = factor*xnorm
1855  if (delta == 0.0d+00) then
1856  delta = factor
1857  end if
1858 
1859  end if
1860 !
1861 ! Form Q'*FVEC and store in QTF.
1862 !
1863  qtf(1:n) = fvec(1:n)
1864 
1865  do j = 1, n
1866  if (fjac(j, j) /= 0.0d+00) then
1867  sum2 = 0.0d+00
1868  do i = j, n
1869  sum2 = sum2 + fjac(i, j)*qtf(i)
1870  end do
1871  temp = -sum2/fjac(j, j)
1872  do i = j, n
1873  qtf(i) = qtf(i) + fjac(i, j)*temp
1874  end do
1875  end if
1876  end do
1877 !
1878 ! Copy the triangular factor of the QR factorization into R.
1879 !
1880  sing = .false.
1881 
1882  do j = 1, n
1883  l = j
1884  do i = 1, j - 1
1885  r(l) = fjac(i, j)
1886  l = l + n - i
1887  end do
1888  r(l) = wa1(j)
1889  if (wa1(j) == 0.0d+00) then
1890  sing = .true.
1891  end if
1892  end do
1893 !
1894 ! Accumulate the orthogonal factor in FJAC.
1895 !
1896  call qform(n, n, fjac, ldfjac)
1897 !
1898 ! Rescale if necessary.
1899 !
1900  if (mode /= 2) then
1901  do j = 1, n
1902  diag(j) = max(diag(j), wa2(j))
1903  end do
1904  end if
1905 !
1906 ! Beginning of the inner loop.
1907 !
1908  do
1909 !
1910 ! If requested, call FCN to enable printing of iterates.
1911 !
1912  if (0 < nprint) then
1913 
1914  iflag = 0
1915  if (mod(iter - 1, nprint) == 0) then
1916  call fcn(n, x, fvec, fjac, ldfjac, iflag)
1917  end if
1918 
1919  if (iflag < 0) then
1920  info = iflag
1921  iflag = 0
1922  if (0 < nprint) then
1923  call fcn(n, x, fvec, fjac, ldfjac, iflag)
1924  end if
1925  return
1926  end if
1927 
1928  end if
1929 !
1930 ! Determine the direction P.
1931 !
1932  call dogleg(n, r, lr, diag, qtf, delta, wa1)
1933 !
1934 ! Store the direction P and X + P.
1935 ! Calculate the norm of P.
1936 !
1937  wa1(1:n) = -wa1(1:n)
1938  wa2(1:n) = x(1:n) + wa1(1:n)
1939  wa3(1:n) = diag(1:n)*wa1(1:n)
1940  pnorm = enorm(n, wa3)
1941 !
1942 ! On the first iteration, adjust the initial step bound.
1943 !
1944  if (iter == 1) then
1945  delta = min(delta, pnorm)
1946  end if
1947 !
1948 ! Evaluate the function at X + P and calculate its norm.
1949 !
1950  iflag = 1
1951  call fcn(n, wa2, wa4, fjac, ldfjac, iflag)
1952  nfev = nfev + 1
1953 
1954  if (iflag < 0) then
1955  info = iflag
1956  iflag = 0
1957  if (0 < nprint) then
1958  call fcn(n, x, fvec, fjac, ldfjac, iflag)
1959  end if
1960  return
1961  end if
1962 
1963  fnorm1 = enorm(n, wa4)
1964 !
1965 ! Compute the scaled actual reduction.
1966 !
1967  actred = -1.0d+00
1968  if (fnorm1 < fnorm) then
1969  actred = 1.0d+00 - (fnorm1/fnorm)**2
1970  end if
1971 !
1972 ! Compute the scaled predicted reduction.
1973 !
1974  l = 1
1975  do i = 1, n
1976  sum2 = 0.0d+00
1977  do j = i, n
1978  sum2 = sum2 + r(l)*wa1(j)
1979  l = l + 1
1980  end do
1981  wa3(i) = qtf(i) + sum2
1982  end do
1983 
1984  temp = enorm(n, wa3)
1985  prered = 0.0d+00
1986  if (temp < fnorm) then
1987  prered = 1.0d+00 - (temp/fnorm)**2
1988  end if
1989 !
1990 ! Compute the ratio of the actual to the predicted reduction.
1991 !
1992  if (0.0d+00 < prered) then
1993  ratio = actred/prered
1994  else
1995  ratio = 0.0d+00
1996  end if
1997 !
1998 ! Update the step bound.
1999 !
2000  if (ratio < 0.1d+00) then
2001 
2002  ncsuc = 0
2003  ncfail = ncfail + 1
2004  delta = 0.5d+00*delta
2005 
2006  else
2007 
2008  ncfail = 0
2009  ncsuc = ncsuc + 1
2010 
2011  if (0.5d+00 <= ratio .or. 1 < ncsuc) then
2012  delta = max(delta, pnorm/0.5d+00)
2013  end if
2014 
2015  if (abs(ratio - 1.0d+00) <= 0.1d+00) then
2016  delta = pnorm/0.5d+00
2017  end if
2018 
2019  end if
2020 !
2021 ! Test for successful iteration.
2022 !
2023 
2024 !
2025 ! Successful iteration.
2026 ! Update X, FVEC, and their norms.
2027 !
2028  if (0.0001d+00 <= ratio) then
2029  x(1:n) = wa2(1:n)
2030  wa2(1:n) = diag(1:n)*x(1:n)
2031  fvec(1:n) = wa4(1:n)
2032  xnorm = enorm(n, wa2)
2033  fnorm = fnorm1
2034  iter = iter + 1
2035  end if
2036 !
2037 ! Determine the progress of the iteration.
2038 !
2039  nslow1 = nslow1 + 1
2040  if (0.001d+00 <= actred) then
2041  nslow1 = 0
2042  end if
2043 
2044  if (jeval) then
2045  nslow2 = nslow2 + 1
2046  end if
2047 
2048  if (0.1d+00 <= actred) then
2049  nslow2 = 0
2050  end if
2051 !
2052 ! Test for convergence.
2053 !
2054  if (delta <= xtol*xnorm .or. fnorm == 0.0d+00) then
2055  info = 1
2056  end if
2057 
2058  if (info /= 0) then
2059  iflag = 0
2060  if (0 < nprint) then
2061  call fcn(n, x, fvec, fjac, ldfjac, iflag)
2062  end if
2063  return
2064  end if
2065 !
2066 ! Tests for termination and stringent tolerances.
2067 !
2068  if (maxfev <= nfev) then
2069  info = 2
2070  end if
2071 
2072  if (0.1d+00*max(0.1d+00*delta, pnorm) <= epsmch*xnorm) then
2073  info = 3
2074  end if
2075 
2076  if (nslow2 == 5) then
2077  info = 4
2078  end if
2079 
2080  if (nslow1 == 10) then
2081  info = 5
2082  end if
2083 
2084  if (info /= 0) then
2085  iflag = 0
2086  if (0 < nprint) then
2087  call fcn(n, x, fvec, fjac, ldfjac, iflag)
2088  end if
2089  return
2090  end if
2091 !
2092 ! Criterion for recalculating jacobian.
2093 !
2094  if (ncfail == 2) then
2095  exit
2096  end if
2097 !
2098 ! Calculate the rank one modification to the jacobian
2099 ! and update QTF if necessary.
2100 !
2101  do j = 1, n
2102  sum2 = dot_product(wa4(1:n), fjac(1:n, j))
2103  wa2(j) = (sum2 - wa3(j))/pnorm
2104  wa1(j) = diag(j)*((diag(j)*wa1(j))/pnorm)
2105  if (0.0001d+00 <= ratio) then
2106  qtf(j) = sum2
2107  end if
2108  end do
2109 !
2110 ! Compute the QR factorization of the updated jacobian.
2111 !
2112  call r1updt(n, n, r, lr, wa1, wa2, wa3, sing)
2113  call r1mpyq(n, n, fjac, ldfjac, wa2, wa3)
2114  call r1mpyq(1, n, qtf, 1, wa2, wa3)
2115 !
2116 ! End of the inner loop.
2117 !
2118  jeval = .false.
2119 
2120  end do
2121 !
2122 ! End of the outer loop.
2123 !
2124  end do
2125 
subroutine dogleg(n, r, lr, diag, qtb, delta, x)
subroutine r1updt(m, n, s, ls, u, v, w, sing)
real(kind=8) function enorm(n, x)
subroutine qrfac(m, n, a, lda, pivot, ipvt, lipvt, rdiag, acnorm)
subroutine r1mpyq(m, n, a, lda, v, w)
subroutine qform(m, n, q, ldq)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ hybrj1()

subroutine hybrj1 ( external  fcn,
integer(kind=4)  n,
real(kind=8), dimension(n)  x,
real(kind=8), dimension(n)  fvec,
real(kind=8), dimension(ldfjac, n)  fjac,
integer(kind=4)  ldfjac,
real(kind=8)  tol,
integer(kind=4)  info 
)

Definition at line 2128 of file suews_util_minpack.f95.

References hybrj().

2128 
2129 !*****************************************************************************80
2130 !
2131 !! HYBRJ1 seeks a zero of N equations in N variables by Powell's method.
2132 !
2133 ! Discussion:
2134 !
2135 ! HYBRJ1 finds a zero of a system of N nonlinear functions in N variables
2136 ! by a modification of the Powell hybrid method. This is done by using the
2137 ! more general nonlinear equation solver HYBRJ. The user
2138 ! must provide a subroutine which calculates the functions
2139 ! and the jacobian.
2140 !
2141 ! Licensing:
2142 !
2143 ! This code is distributed under the GNU LGPL license.
2144 !
2145 ! Modified:
2146 !
2147 ! 06 April 2010
2148 !
2149 ! Author:
2150 !
2151 ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom.
2152 ! FORTRAN90 version by John Burkardt.
2153 !
2154 ! Reference:
2155 !
2156 ! Jorge More, Burton Garbow, Kenneth Hillstrom,
2157 ! User Guide for MINPACK-1,
2158 ! Technical Report ANL-80-74,
2159 ! Argonne National Laboratory, 1980.
2160 !
2161 ! Parameters:
2162 !
2163 ! Input, external FCN, the name of the user-supplied subroutine which
2164 ! calculates the functions and the jacobian. FCN should have the form:
2165 ! subroutine fcn ( n, x, fvec, fjac, ldfjac, iflag )
2166 ! integer ( kind = 4 ) ldfjac
2167 ! integer ( kind = 4 ) n
2168 ! real ( kind = 8 ) fjac(ldfjac,n)
2169 ! real ( kind = 8 ) fvec(n)
2170 ! integer ( kind = 4 ) iflag
2171 ! real ( kind = 8 ) x(n)
2172 !
2173 ! If IFLAG = 0 on input, then FCN is only being called to allow the user
2174 ! to print out the current iterate.
2175 ! If IFLAG = 1 on input, FCN should calculate the functions at X and
2176 ! return this vector in FVEC.
2177 ! If IFLAG = 2 on input, FCN should calculate the jacobian at X and
2178 ! return this matrix in FJAC.
2179 ! To terminate the algorithm, FCN may set IFLAG negative on return.
2180 !
2181 ! Input, integer ( kind = 4 ) N, the number of functions and variables.
2182 !
2183 ! Input/output, real ( kind = 8 ) X(N). On input, X must contain an initial
2184 ! estimate of the solution vector. On output X contains the final
2185 ! estimate of the solution vector.
2186 !
2187 ! Output, real ( kind = 8 ) FVEC(N), the functions evaluated at the output X.
2188 !
2189 ! Output, real ( kind = 8 ) FJAC(LDFJAC,N), an N by N array which contains
2190 ! the orthogonal matrix Q produced by the QR factorization of the final
2191 ! approximate jacobian.
2192 !
2193 ! Input, integer ( kind = 4 ) LDFJAC, the leading dimension of FJAC.
2194 ! LDFJAC must be at least N.
2195 !
2196 ! Input, real ( kind = 8 ) TOL. Termination occurs when the algorithm
2197 ! estimates that the relative error between X and the solution is at most
2198 ! TOL. TOL should be nonnegative.
2199 !
2200 ! Output, integer ( kind = 4 ) INFO, error flag. If the user has terminated
2201 ! execution, INFO is set to the (negative) value of IFLAG. See description
2202 ! of FCN. Otherwise, INFO is set as follows:
2203 ! 0, improper input parameters.
2204 ! 1, algorithm estimates that the relative error between X and the
2205 ! solution is at most TOL.
2206 ! 2, number of calls to FCN with IFLAG = 1 has reached 100*(N+1).
2207 ! 3, TOL is too small. No further improvement in the approximate
2208 ! solution X is possible.
2209 ! 4, iteration is not making good progress.
2210 !
2211  implicit none
2212 
2213  integer(kind=4) ldfjac
2214  integer(kind=4) n
2215 
2216  real(kind=8) diag(n)
2217  real(kind=8) factor
2218  external fcn
2219  real(kind=8) fjac(ldfjac, n)
2220  real(kind=8) fvec(n)
2221  integer(kind=4) info
2222  ! integer ( kind = 4 ) j
2223  integer(kind=4) lr
2224  integer(kind=4) maxfev
2225  integer(kind=4) mode
2226  integer(kind=4) nfev
2227  integer(kind=4) njev
2228  integer(kind=4) nprint
2229  real(kind=8) qtf(n)
2230  real(kind=8) r((n*(n + 1))/2)
2231  real(kind=8) tol
2232  real(kind=8) x(n)
2233  real(kind=8) xtol
2234 
2235  info = 0
2236 
2237  if (n <= 0) then
2238  return
2239  else if (ldfjac < n) then
2240  return
2241  else if (tol < 0.0d+00) then
2242  return
2243  end if
2244 
2245  maxfev = 100*(n + 1)
2246  xtol = tol
2247  mode = 2
2248  diag(1:n) = 1.0d+00
2249  factor = 100.0d+00
2250  nprint = 0
2251  lr = (n*(n + 1))/2
2252 
2253  call hybrj(fcn, n, x, fvec, fjac, ldfjac, xtol, maxfev, diag, mode, &
2254  factor, nprint, info, nfev, njev, r, lr, qtf)
2255 
2256  if (info == 5) then
2257  info = 4
2258  end if
2259 
2260  return
subroutine hybrj(fcn, n, x, fvec, fjac, ldfjac, xtol, maxfev, diag, mode, factor, nprint, info, nfev, njev, r, lr, qtf)
Here is the call graph for this function:

◆ lmder()

subroutine lmder ( external  fcn,
integer(kind=4)  m,
integer(kind=4)  n,
real(kind=8), dimension(n)  x,
real(kind=8), dimension(m)  fvec,
real(kind=8), dimension(ldfjac, n)  fjac,
integer(kind=4)  ldfjac,
real(kind=8)  ftol,
real(kind=8)  xtol,
real(kind=8)  gtol,
integer(kind=4)  maxfev,
real(kind=8), dimension(n)  diag,
integer(kind=4)  mode,
real(kind=8)  factor,
integer(kind=4)  nprint,
integer(kind=4)  info,
integer(kind=4)  nfev,
integer(kind=4)  njev,
integer(kind=4), dimension(n)  ipvt,
real(kind=8), dimension(n)  qtf 
)

Definition at line 2264 of file suews_util_minpack.f95.

References lmpar(), and qrfac().

Referenced by lmder1().

2264 
2265 !*****************************************************************************80
2266 !
2267 !! LMDER minimizes M functions in N variables by the Levenberg-Marquardt method.
2268 !
2269 ! Discussion:
2270 !
2271 ! LMDER minimizes the sum of the squares of M nonlinear functions in
2272 ! N variables by a modification of the Levenberg-Marquardt algorithm.
2273 ! The user must provide a subroutine which calculates the functions
2274 ! and the jacobian.
2275 !
2276 ! Licensing:
2277 !
2278 ! This code is distributed under the GNU LGPL license.
2279 !
2280 ! Modified:
2281 !
2282 ! 06 April 2010
2283 !
2284 ! Author:
2285 !
2286 ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom.
2287 ! FORTRAN90 version by John Burkardt.
2288 !
2289 ! Reference:
2290 !
2291 ! Jorge More, Burton Garbow, Kenneth Hillstrom,
2292 ! User Guide for MINPACK-1,
2293 ! Technical Report ANL-80-74,
2294 ! Argonne National Laboratory, 1980.
2295 !
2296 ! Parameters:
2297 !
2298 ! Input, external FCN, the name of the user-supplied subroutine which
2299 ! calculates the functions and the jacobian. FCN should have the form:
2300 ! subroutine fcn ( m, n, x, fvec, fjac, ldfjac, iflag )
2301 ! integer ( kind = 4 ) ldfjac
2302 ! integer ( kind = 4 ) n
2303 ! real ( kind = 8 ) fjac(ldfjac,n)
2304 ! real ( kind = 8 ) fvec(m)
2305 ! integer ( kind = 4 ) iflag
2306 ! real ( kind = 8 ) x(n)
2307 !
2308 ! If IFLAG = 0 on input, then FCN is only being called to allow the user
2309 ! to print out the current iterate.
2310 ! If IFLAG = 1 on input, FCN should calculate the functions at X and
2311 ! return this vector in FVEC.
2312 ! If IFLAG = 2 on input, FCN should calculate the jacobian at X and
2313 ! return this matrix in FJAC.
2314 ! To terminate the algorithm, FCN may set IFLAG negative on return.
2315 !
2316 ! Input, integer ( kind = 4 ) M, is the number of functions.
2317 !
2318 ! Input, integer ( kind = 4 ) N, is the number of variables.
2319 ! N must not exceed M.
2320 !
2321 ! Input/output, real ( kind = 8 ) X(N). On input, X must contain an initial
2322 ! estimate of the solution vector. On output X contains the final
2323 ! estimate of the solution vector.
2324 !
2325 ! Output, real ( kind = 8 ) FVEC(M), the functions evaluated at the output X.
2326 !
2327 ! Output, real ( kind = 8 ) FJAC(LDFJAC,N), an M by N array. The upper
2328 ! N by N submatrix of FJAC contains an upper triangular matrix R with
2329 ! diagonal elements of nonincreasing magnitude such that
2330 ! P' * ( JAC' * JAC ) * P = R' * R,
2331 ! where P is a permutation matrix and JAC is the final calculated jacobian.
2332 ! Column J of P is column IPVT(J) of the identity matrix. The lower
2333 ! trapezoidal part of FJAC contains information generated during
2334 ! the computation of R.
2335 !
2336 ! Input, integer ( kind = 4 ) LDFJAC, the leading dimension of FJAC.
2337 ! LDFJAC must be at least M.
2338 !
2339 ! Input, real ( kind = 8 ) FTOL. Termination occurs when both the actual
2340 ! and predicted relative reductions in the sum of squares are at most FTOL.
2341 ! Therefore, FTOL measures the relative error desired in the sum of
2342 ! squares. FTOL should be nonnegative.
2343 !
2344 ! Input, real ( kind = 8 ) XTOL. Termination occurs when the relative error
2345 ! between two consecutive iterates is at most XTOL. XTOL should be
2346 ! nonnegative.
2347 !
2348 ! Input, real ( kind = 8 ) GTOL. Termination occurs when the cosine of the
2349 ! angle between FVEC and any column of the jacobian is at most GTOL in
2350 ! absolute value. Therefore, GTOL measures the orthogonality desired
2351 ! between the function vector and the columns of the jacobian. GTOL should
2352 ! be nonnegative.
2353 !
2354 ! Input, integer ( kind = 4 ) MAXFEV. Termination occurs when the number of
2355 ! calls to FCN with IFLAG = 1 is at least MAXFEV by the end of an iteration.
2356 !
2357 ! Input/output, real ( kind = 8 ) DIAG(N). If MODE = 1, then DIAG is set
2358 ! internally. If MODE = 2, then DIAG must contain positive entries that
2359 ! serve as multiplicative scale factors for the variables.
2360 !
2361 ! Input, integer ( kind = 4 ) MODE, scaling option.
2362 ! 1, variables will be scaled internally.
2363 ! 2, scaling is specified by the input DIAG vector.
2364 !
2365 ! Input, real ( kind = 8 ) FACTOR, determines the initial step bound. This
2366 ! bound is set to the product of FACTOR and the euclidean norm of DIAG*X if
2367 ! nonzero, or else to FACTOR itself. In most cases, FACTOR should lie
2368 ! in the interval (0.1, 100) with 100 the recommended value.
2369 !
2370 ! Input, integer ( kind = 4 ) NPRINT, enables controlled printing of iterates
2371 ! if it is positive. In this case, FCN is called with IFLAG = 0 at the
2372 ! beginning of the first iteration and every NPRINT iterations thereafter
2373 ! and immediately prior to return, with X and FVEC available
2374 ! for printing. If NPRINT is not positive, no special calls
2375 ! of FCN with IFLAG = 0 are made.
2376 !
2377 ! Output, integer ( kind = 4 ) INFO, error flag. If the user has terminated
2378 ! execution, INFO is set to the (negative) value of IFLAG. See description
2379 ! of FCN. Otherwise, INFO is set as follows:
2380 ! 0, improper input parameters.
2381 ! 1, both actual and predicted relative reductions in the sum of
2382 ! squares are at most FTOL.
2383 ! 2, relative error between two consecutive iterates is at most XTOL.
2384 ! 3, conditions for INFO = 1 and INFO = 2 both hold.
2385 ! 4, the cosine of the angle between FVEC and any column of the jacobian
2386 ! is at most GTOL in absolute value.
2387 ! 5, number of calls to FCN with IFLAG = 1 has reached MAXFEV.
2388 ! 6, FTOL is too small. No further reduction in the sum of squares
2389 ! is possible.
2390 ! 7, XTOL is too small. No further improvement in the approximate
2391 ! solution X is possible.
2392 ! 8, GTOL is too small. FVEC is orthogonal to the columns of the
2393 ! jacobian to machine precision.
2394 !
2395 ! Output, integer ( kind = 4 ) NFEV, the number of calls to FCN with
2396 ! IFLAG = 1.
2397 !
2398 ! Output, integer ( kind = 4 ) NJEV, the number of calls to FCN with
2399 ! IFLAG = 2.
2400 !
2401 ! Output, integer ( kind = 4 ) IPVT(N), defines a permutation matrix P
2402 ! such that JAC*P = Q*R, where JAC is the final calculated jacobian, Q is
2403 ! orthogonal (not stored), and R is upper triangular with diagonal
2404 ! elements of nonincreasing magnitude. Column J of P is column
2405 ! IPVT(J) of the identity matrix.
2406 !
2407 ! Output, real ( kind = 8 ) QTF(N), contains the first N elements of Q'*FVEC.
2408 !
2409  implicit none
2410 
2411  integer(kind=4) ldfjac
2412  integer(kind=4) m
2413  integer(kind=4) n
2414 
2415  real(kind=8) actred
2416  real(kind=8) delta
2417  real(kind=8) diag(n)
2418  real(kind=8) dirder
2419  real(kind=8) enorm
2420  real(kind=8) epsmch
2421  real(kind=8) factor
2422  external fcn
2423  real(kind=8) fjac(ldfjac, n)
2424  real(kind=8) fnorm
2425  real(kind=8) fnorm1
2426  real(kind=8) ftol
2427  real(kind=8) fvec(m)
2428  real(kind=8) gnorm
2429  real(kind=8) gtol
2430  ! integer ( kind = 4 ) i
2431  integer(kind=4) iflag
2432  integer(kind=4) info
2433  integer(kind=4) ipvt(n)
2434  integer(kind=4) iter
2435  integer(kind=4) j
2436  integer(kind=4) l
2437  integer(kind=4) maxfev
2438  integer(kind=4) mode
2439  integer(kind=4) nfev
2440  integer(kind=4) njev
2441  integer(kind=4) nprint
2442  real(kind=8) par
2443  logical pivot
2444  real(kind=8) pnorm
2445  real(kind=8) prered
2446  real(kind=8) qtf(n)
2447  real(kind=8) ratio
2448  real(kind=8) sum2
2449  real(kind=8) temp
2450  real(kind=8) temp1
2451  real(kind=8) temp2
2452  real(kind=8) wa1(n)
2453  real(kind=8) wa2(n)
2454  real(kind=8) wa3(n)
2455  real(kind=8) wa4(m)
2456  real(kind=8) xnorm
2457  real(kind=8) x(n)
2458  real(kind=8) xtol
2459 
2460  epsmch = epsilon(epsmch)
2461 
2462  info = 0
2463  iflag = 0
2464  nfev = 0
2465  njev = 0
2466 !
2467 ! Check the input parameters for errors.
2468 !
2469  if (n <= 0) then
2470  go to 300
2471  end if
2472 
2473  if (m < n) then
2474  go to 300
2475  end if
2476 
2477  if (ldfjac < m &
2478  .or. ftol < 0.0d+00 .or. xtol < 0.0d+00 .or. gtol < 0.0d+00 &
2479  .or. maxfev <= 0 .or. factor <= 0.0d+00) then
2480  go to 300
2481  end if
2482 
2483  if (mode == 2) then
2484  do j = 1, n
2485  if (diag(j) <= 0.0d+00) then
2486  go to 300
2487  end if
2488  end do
2489  end if
2490 !
2491 ! Evaluate the function at the starting point and calculate its norm.
2492 !
2493  iflag = 1
2494  call fcn(m, n, x, fvec, fjac, ldfjac, iflag)
2495  nfev = 1
2496  if (iflag < 0) then
2497  go to 300
2498  end if
2499 
2500  fnorm = enorm(m, fvec)
2501 !
2502 ! Initialize Levenberg-Marquardt parameter and iteration counter.
2503 !
2504  par = 0.0d+00
2505  iter = 1
2506 !
2507 ! Beginning of the outer loop.
2508 !
2509 30 continue
2510 !
2511 ! Calculate the jacobian matrix.
2512 !
2513  iflag = 2
2514  call fcn(m, n, x, fvec, fjac, ldfjac, iflag)
2515 
2516  njev = njev + 1
2517 
2518  if (iflag < 0) then
2519  go to 300
2520  end if
2521 !
2522 ! If requested, call FCN to enable printing of iterates.
2523 !
2524  if (0 < nprint) then
2525  iflag = 0
2526  if (mod(iter - 1, nprint) == 0) then
2527  call fcn(m, n, x, fvec, fjac, ldfjac, iflag)
2528  end if
2529  if (iflag < 0) then
2530  go to 300
2531  end if
2532  end if
2533 !
2534 ! Compute the QR factorization of the jacobian.
2535 !
2536  pivot = .true.
2537  call qrfac(m, n, fjac, ldfjac, pivot, ipvt, n, wa1, wa2)
2538 !
2539 ! On the first iteration and if mode is 1, scale according
2540 ! to the norms of the columns of the initial jacobian.
2541 !
2542  if (iter == 1) then
2543 
2544  if (mode /= 2) then
2545  diag(1:n) = wa2(1:n)
2546  do j = 1, n
2547  if (wa2(j) == 0.0d+00) then
2548  diag(j) = 1.0d+00
2549  end if
2550  end do
2551  end if
2552 !
2553 ! On the first iteration, calculate the norm of the scaled X
2554 ! and initialize the step bound DELTA.
2555 !
2556  wa3(1:n) = diag(1:n)*x(1:n)
2557 
2558  xnorm = enorm(n, wa3)
2559  delta = factor*xnorm
2560  if (delta == 0.0d+00) then
2561  delta = factor
2562  end if
2563 
2564  end if
2565 !
2566 ! Form Q'*FVEC and store the first N components in QTF.
2567 !
2568  wa4(1:m) = fvec(1:m)
2569 
2570  do j = 1, n
2571 
2572  if (fjac(j, j) /= 0.0d+00) then
2573  sum2 = dot_product(wa4(j:m), fjac(j:m, j))
2574  temp = -sum2/fjac(j, j)
2575  wa4(j:m) = wa4(j:m) + fjac(j:m, j)*temp
2576  end if
2577 
2578  fjac(j, j) = wa1(j)
2579  qtf(j) = wa4(j)
2580 
2581  end do
2582 !
2583 ! Compute the norm of the scaled gradient.
2584 !
2585  gnorm = 0.0d+00
2586 
2587  if (fnorm /= 0.0d+00) then
2588 
2589  do j = 1, n
2590  l = ipvt(j)
2591  if (wa2(l) /= 0.0d+00) then
2592  sum2 = dot_product(qtf(1:j), fjac(1:j, j))/fnorm
2593  gnorm = max(gnorm, abs(sum2/wa2(l)))
2594  end if
2595  end do
2596 
2597  end if
2598 !
2599 ! Test for convergence of the gradient norm.
2600 !
2601  if (gnorm <= gtol) then
2602  info = 4
2603  go to 300
2604  end if
2605 !
2606 ! Rescale if necessary.
2607 !
2608  if (mode /= 2) then
2609  do j = 1, n
2610  diag(j) = max(diag(j), wa2(j))
2611  end do
2612  end if
2613 !
2614 ! Beginning of the inner loop.
2615 !
2616 200 continue
2617 !
2618 ! Determine the Levenberg-Marquardt parameter.
2619 !
2620  call lmpar(n, fjac, ldfjac, ipvt, diag, qtf, delta, par, wa1, wa2)
2621 !
2622 ! Store the direction p and x + p. calculate the norm of p.
2623 !
2624  wa1(1:n) = -wa1(1:n)
2625  wa2(1:n) = x(1:n) + wa1(1:n)
2626  wa3(1:n) = diag(1:n)*wa1(1:n)
2627 
2628  pnorm = enorm(n, wa3)
2629 !
2630 ! On the first iteration, adjust the initial step bound.
2631 !
2632  if (iter == 1) then
2633  delta = min(delta, pnorm)
2634  end if
2635 !
2636 ! Evaluate the function at x + p and calculate its norm.
2637 !
2638  iflag = 1
2639  call fcn(m, n, wa2, wa4, fjac, ldfjac, iflag)
2640 
2641  nfev = nfev + 1
2642 
2643  if (iflag < 0) then
2644  go to 300
2645  end if
2646 
2647  fnorm1 = enorm(m, wa4)
2648 !
2649 ! Compute the scaled actual reduction.
2650 !
2651  actred = -1.0d+00
2652  if (0.1d+00*fnorm1 < fnorm) then
2653  actred = 1.0d+00 - (fnorm1/fnorm)**2
2654  end if
2655 !
2656 ! Compute the scaled predicted reduction and
2657 ! the scaled directional derivative.
2658 !
2659  do j = 1, n
2660  wa3(j) = 0.0d+00
2661  l = ipvt(j)
2662  temp = wa1(l)
2663  wa3(1:j) = wa3(1:j) + fjac(1:j, j)*temp
2664  end do
2665 
2666  temp1 = enorm(n, wa3)/fnorm
2667  temp2 = (sqrt(par)*pnorm)/fnorm
2668  prered = temp1**2 + temp2**2/0.5d+00
2669  dirder = -(temp1**2 + temp2**2)
2670 !
2671 ! Compute the ratio of the actual to the predicted reduction.
2672 !
2673  if (prered /= 0.0d+00) then
2674  ratio = actred/prered
2675  else
2676  ratio = 0.0d+00
2677  end if
2678 !
2679 ! Update the step bound.
2680 !
2681  if (ratio <= 0.25d+00) then
2682 
2683  if (0.0d+00 <= actred) then
2684  temp = 0.5d+00
2685  end if
2686 
2687  if (actred < 0.0d+00) then
2688  temp = 0.5d+00*dirder/(dirder + 0.5d+00*actred)
2689  end if
2690 
2691  if (0.1d+00*fnorm1 >= fnorm .or. temp < 0.1d+00) then
2692  temp = 0.1d+00
2693  end if
2694 
2695  delta = temp*min(delta, pnorm/0.1d+00)
2696  par = par/temp
2697 
2698  else
2699 
2700  if (par == 0.0d+00 .or. ratio >= 0.75d+00) then
2701  delta = 2.0d+00*pnorm
2702  par = 0.5d+00*par
2703  end if
2704 
2705  end if
2706 !
2707 ! Successful iteration.
2708 !
2709 ! Update X, FVEC, and their norms.
2710 !
2711  if (0.0001d+00 <= ratio) then
2712  x(1:n) = wa2(1:n)
2713  wa2(1:n) = diag(1:n)*x(1:n)
2714  fvec(1:m) = wa4(1:m)
2715  xnorm = enorm(n, wa2)
2716  fnorm = fnorm1
2717  iter = iter + 1
2718  end if
2719 !
2720 ! Tests for convergence.
2721 !
2722  if (abs(actred) <= ftol .and. &
2723  prered <= ftol .and. &
2724  0.5d+00*ratio <= 1.0d+00) then
2725  info = 1
2726  end if
2727 
2728  if (delta <= xtol*xnorm) then
2729  info = 2
2730  end if
2731 
2732  if (abs(actred) <= ftol .and. prered <= ftol &
2733  .and. 0.5d+00*ratio <= 1.0d+00 .and. info == 2) then
2734  info = 3
2735  end if
2736 
2737  if (info /= 0) then
2738  go to 300
2739  end if
2740 !
2741 ! Tests for termination and stringent tolerances.
2742 !
2743  if (nfev >= maxfev) then
2744  info = 5
2745  end if
2746 
2747  if (abs(actred) <= epsmch .and. prered <= epsmch &
2748  .and. 0.5d+00*ratio <= 1.0d+00) then
2749  info = 6
2750  end if
2751 
2752  if (delta <= epsmch*xnorm) then
2753  info = 7
2754  end if
2755 
2756  if (gnorm <= epsmch) then
2757  info = 8
2758  end if
2759 
2760  if (info /= 0) then
2761  go to 300
2762  end if
2763 !
2764 ! End of the inner loop. repeat if iteration unsuccessful.
2765 !
2766  if (ratio < 0.0001d+00) then
2767  go to 200
2768  end if
2769 !
2770 ! End of the outer loop.
2771 !
2772  go to 30
2773 
2774 300 continue
2775 !
2776 ! Termination, either normal or user imposed.
2777 !
2778  if (iflag < 0) then
2779  info = iflag
2780  end if
2781 
2782  iflag = 0
2783 
2784  if (0 < nprint) then
2785  call fcn(m, n, x, fvec, fjac, ldfjac, iflag)
2786  end if
2787 
2788  return
real(kind=8) function enorm(n, x)
subroutine qrfac(m, n, a, lda, pivot, ipvt, lipvt, rdiag, acnorm)
subroutine lmpar(n, r, ldr, ipvt, diag, qtb, delta, par, x, sdiag)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ lmder1()

subroutine lmder1 ( external  fcn,
integer(kind=4)  m,
integer(kind=4)  n,
real(kind=8), dimension(n)  x,
real(kind=8), dimension(m)  fvec,
real(kind=8), dimension(ldfjac, n)  fjac,
integer(kind=4)  ldfjac,
real(kind=8)  tol,
integer(kind=4)  info 
)

Definition at line 2791 of file suews_util_minpack.f95.

References lmder().

2791 
2792 !*****************************************************************************80
2793 !
2794 !! LMDER1 minimizes M functions in N variables by Levenberg-Marquardt method.
2795 !
2796 ! Discussion:
2797 !
2798 ! LMDER1 minimizes the sum of the squares of M nonlinear functions in
2799 ! N variables by a modification of the Levenberg-Marquardt algorithm.
2800 ! This is done by using the more general least-squares solver LMDER.
2801 ! The user must provide a subroutine which calculates the functions
2802 ! and the jacobian.
2803 !
2804 ! Licensing:
2805 !
2806 ! This code is distributed under the GNU LGPL license.
2807 !
2808 ! Modified:
2809 !
2810 ! 06 April 2010
2811 !
2812 ! Author:
2813 !
2814 ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom.
2815 ! FORTRAN90 version by John Burkardt.
2816 !
2817 ! Reference:
2818 !
2819 ! Jorge More, Burton Garbow, Kenneth Hillstrom,
2820 ! User Guide for MINPACK-1,
2821 ! Technical Report ANL-80-74,
2822 ! Argonne National Laboratory, 1980.
2823 !
2824 ! Parameters:
2825 !
2826 ! Input, external FCN, the name of the user-supplied subroutine which
2827 ! calculates the functions and the jacobian. FCN should have the form:
2828 ! subroutine fcn ( m, n, x, fvec, fjac, ldfjac, iflag )
2829 ! integer ( kind = 4 ) ldfjac
2830 ! integer ( kind = 4 ) n
2831 ! real ( kind = 8 ) fjac(ldfjac,n)
2832 ! real ( kind = 8 ) fvec(m)
2833 ! integer ( kind = 4 ) iflag
2834 ! real ( kind = 8 ) x(n)
2835 !
2836 ! If IFLAG = 0 on input, then FCN is only being called to allow the user
2837 ! to print out the current iterate.
2838 ! If IFLAG = 1 on input, FCN should calculate the functions at X and
2839 ! return this vector in FVEC.
2840 ! If IFLAG = 2 on input, FCN should calculate the jacobian at X and
2841 ! return this matrix in FJAC.
2842 ! To terminate the algorithm, FCN may set IFLAG negative on return.
2843 !
2844 ! Input, integer ( kind = 4 ) M, the number of functions.
2845 !
2846 ! Input, integer ( kind = 4 ) N, is the number of variables.
2847 ! N must not exceed M.
2848 !
2849 ! Input/output, real ( kind = 8 ) X(N). On input, X must contain an initial
2850 ! estimate of the solution vector. On output X contains the final
2851 ! estimate of the solution vector.
2852 !
2853 ! Output, real ( kind = 8 ) FVEC(M), the functions evaluated at the output X.
2854 !
2855 ! Output, real ( kind = 8 ) FJAC(LDFJAC,N), an M by N array. The upper
2856 ! N by N submatrix contains an upper triangular matrix R with
2857 ! diagonal elements of nonincreasing magnitude such that
2858 ! P' * ( JAC' * JAC ) * P = R' * R,
2859 ! where P is a permutation matrix and JAC is the final calculated
2860 ! jacobian. Column J of P is column IPVT(J) of the identity matrix.
2861 ! The lower trapezoidal part of FJAC contains information generated during
2862 ! the computation of R.
2863 !
2864 ! Input, integer ( kind = 4 ) LDFJAC, is the leading dimension of FJAC,
2865 ! which must be no less than M.
2866 !
2867 ! Input, real ( kind = 8 ) TOL. Termination occurs when the algorithm
2868 ! estimates either that the relative error in the sum of squares is at
2869 ! most TOL or that the relative error between X and the solution is at
2870 ! most TOL.
2871 !
2872 ! Output, integer ( kind = 4 ) INFO, error flag. If the user has terminated
2873 ! execution, INFO is set to the (negative) value of IFLAG. See description
2874 ! of FCN. Otherwise, INFO is set as follows:
2875 ! 0, improper input parameters.
2876 ! 1, algorithm estimates that the relative error in the sum of squares
2877 ! is at most TOL.
2878 ! 2, algorithm estimates that the relative error between X and the
2879 ! solution is at most TOL.
2880 ! 3, conditions for INFO = 1 and INFO = 2 both hold.
2881 ! 4, FVEC is orthogonal to the columns of the jacobian to machine precision.
2882 ! 5, number of calls to FCN with IFLAG = 1 has reached 100*(N+1).
2883 ! 6, TOL is too small. No further reduction in the sum of squares is
2884 ! possible.
2885 ! 7, TOL is too small. No further improvement in the approximate
2886 ! solution X is possible.
2887 !
2888  implicit none
2889 
2890  integer(kind=4) ldfjac
2891  integer(kind=4) m
2892  integer(kind=4) n
2893 
2894  real(kind=8) diag(n)
2895  real(kind=8) factor
2896  external fcn
2897  real(kind=8) fjac(ldfjac, n)
2898  real(kind=8) ftol
2899  real(kind=8) fvec(m)
2900  real(kind=8) gtol
2901  integer(kind=4) info
2902  integer(kind=4) ipvt(n)
2903  integer(kind=4) maxfev
2904  integer(kind=4) mode
2905  integer(kind=4) nfev
2906  integer(kind=4) njev
2907  integer(kind=4) nprint
2908  real(kind=8) qtf(n)
2909  real(kind=8) tol
2910  real(kind=8) x(n)
2911  real(kind=8) xtol
2912 
2913  info = 0
2914 
2915  if (n <= 0) then
2916  return
2917  else if (m < n) then
2918  return
2919  else if (ldfjac < m) then
2920  return
2921  else if (tol < 0.0d+00) then
2922  return
2923  end if
2924 
2925  factor = 100.0d+00
2926  maxfev = 100*(n + 1)
2927  ftol = tol
2928  xtol = tol
2929  gtol = 0.0d+00
2930  mode = 1
2931  nprint = 0
2932 
2933  call lmder(fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, &
2934  diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf)
2935 
2936  if (info == 8) then
2937  info = 4
2938  end if
2939 
2940  return
subroutine lmder(fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf)
Here is the call graph for this function:

◆ lmdif()

subroutine lmdif ( external  fcn,
integer(kind=4)  m,
integer(kind=4)  n,
real(kind=8), dimension(n)  x,
real(kind=8), dimension(m)  xdat,
real(kind=8), dimension(m)  ydat,
real(kind=8), dimension(m)  fvec,
real(kind=8)  ftol,
real(kind=8)  xtol,
real(kind=8)  gtol,
integer(kind=4)  maxfev,
real(kind=8)  epsfcn,
real(kind=8), dimension(n)  diag,
integer(kind=4)  mode,
real(kind=8)  factor,
integer(kind=4)  nprint,
integer(kind=4)  info,
integer(kind=4)  nfev,
real(kind=8), dimension(ldfjac, n)  fjac,
integer(kind=4)  ldfjac,
integer(kind=4), dimension(n)  ipvt,
real(kind=8), dimension(n)  qtf 
)

Definition at line 2944 of file suews_util_minpack.f95.

References fdjac2(), lmpar(), and qrfac().

Referenced by lmdif1().

2944 
2945 !*****************************************************************************80
2946 !
2947 !! LMDIF minimizes M functions in N variables by the Levenberg-Marquardt method.
2948 !
2949 ! Discussion:
2950 !
2951 ! LMDIF minimizes the sum of the squares of M nonlinear functions in
2952 ! N variables by a modification of the Levenberg-Marquardt algorithm.
2953 ! The user must provide a subroutine which calculates the functions.
2954 ! The jacobian is then calculated by a forward-difference approximation.
2955 !
2956 ! Licensing:
2957 !
2958 ! This code is distributed under the GNU LGPL license.
2959 !
2960 ! Modified:
2961 !
2962 ! 20 Jul 2017
2963 !
2964 ! Author:
2965 !
2966 ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom.
2967 ! FORTRAN90 version by John Burkardt.
2968 ! adapted for AnOHM in SUEWS, Ting Sun, ting.sun@reading.ac.uk
2969 !
2970 ! Reference:
2971 !
2972 ! Jorge More, Burton Garbow, Kenneth Hillstrom,
2973 ! User Guide for MINPACK-1,
2974 ! Technical Report ANL-80-74,
2975 ! Argonne National Laboratory, 1980.
2976 !
2977 ! Parameters:
2978 !
2979 ! Input, external FCN, the name of the user-supplied subroutine which
2980 ! calculates the functions. The routine should have the form:
2981 ! subroutine fcn ( m, n, x, xdat, ydat, fvec, iflag ) ! xdat, ydat added for AnOHM, TS 20 Jul 2017
2982 ! integer ( kind = 4 ) m
2983 ! integer ( kind = 4 ) n
2984 ! real ( kind = 8 ) fvec(m),xdat(m),ydat(m)
2985 ! integer ( kind = 4 ) iflag
2986 ! real ( kind = 8 ) x(n)
2987 ! If IFLAG = 0 on input, then FCN is only being called to allow the user
2988 ! to print out the current iterate.
2989 ! If IFLAG = 1 on input, FCN should calculate the functions at X and
2990 ! return this vector in FVEC.
2991 ! To terminate the algorithm, FCN may set IFLAG negative on return.
2992 !
2993 ! Input, integer ( kind = 4 ) M, the number of functions.
2994 !
2995 ! Input, integer ( kind = 4 ) N, the number of variables.
2996 ! N must not exceed M.
2997 !
2998 ! Input/output, real ( kind = 8 ) X(N). On input, X must contain an initial
2999 ! estimate of the solution vector. On output X contains the final
3000 ! estimate of the solution vector.
3001 
3002 ! Input, real ( kind = 8 ) XDAT(M), YDAT(M). On input, XDAT/YDAT must contain
3003 ! observations to construct FVEC(M). ! added for AnOHM, TS 20 Jul 2017
3004 !
3005 ! Output, real ( kind = 8 ) FVEC(M), the functions evaluated at the output X.
3006 !
3007 ! Input, real ( kind = 8 ) FTOL. Termination occurs when both the actual
3008 ! and predicted relative reductions in the sum of squares are at most FTOL.
3009 ! Therefore, FTOL measures the relative error desired in the sum of
3010 ! squares. FTOL should be nonnegative.
3011 !
3012 ! Input, real ( kind = 8 ) XTOL. Termination occurs when the relative error
3013 ! between two consecutive iterates is at most XTOL. Therefore, XTOL
3014 ! measures the relative error desired in the approximate solution. XTOL
3015 ! should be nonnegative.
3016 !
3017 ! Input, real ( kind = 8 ) GTOL. termination occurs when the cosine of the
3018 ! angle between FVEC and any column of the jacobian is at most GTOL in
3019 ! absolute value. Therefore, GTOL measures the orthogonality desired
3020 ! between the function vector and the columns of the jacobian. GTOL should
3021 ! be nonnegative.
3022 !
3023 ! Input, integer ( kind = 4 ) MAXFEV. Termination occurs when the number of
3024 ! calls to FCN is at least MAXFEV by the end of an iteration.
3025 !
3026 ! Input, real ( kind = 8 ) EPSFCN, is used in determining a suitable step
3027 ! length for the forward-difference approximation. This approximation
3028 ! assumes that the relative errors in the functions are of the order of
3029 ! EPSFCN. If EPSFCN is less than the machine precision, it is assumed that
3030 ! the relative errors in the functions are of the order of the machine
3031 ! precision.
3032 !
3033 ! Input/output, real ( kind = 8 ) DIAG(N). If MODE = 1, then DIAG is set
3034 ! internally. If MODE = 2, then DIAG must contain positive entries that
3035 ! serve as multiplicative scale factors for the variables.
3036 !
3037 ! Input, integer ( kind = 4 ) MODE, scaling option.
3038 ! 1, variables will be scaled internally.
3039 ! 2, scaling is specified by the input DIAG vector.
3040 !
3041 ! Input, real ( kind = 8 ) FACTOR, determines the initial step bound.
3042 ! This bound is set to the product of FACTOR and the euclidean norm of
3043 ! DIAG*X if nonzero, or else to FACTOR itself. In most cases, FACTOR should
3044 ! lie in the interval (0.1, 100) with 100 the recommended value.
3045 !
3046 ! Input, integer ( kind = 4 ) NPRINT, enables controlled printing of iterates
3047 ! if it is positive. In this case, FCN is called with IFLAG = 0 at the
3048 ! beginning of the first iteration and every NPRINT iterations thereafter
3049 ! and immediately prior to return, with X and FVEC available
3050 ! for printing. If NPRINT is not positive, no special calls
3051 ! of FCN with IFLAG = 0 are made.
3052 !
3053 ! Output, integer ( kind = 4 ) INFO, error flag. If the user has terminated
3054 ! execution, INFO is set to the (negative) value of IFLAG. See description
3055 ! of FCN. Otherwise, INFO is set as follows:
3056 ! 0, improper input parameters.
3057 ! 1, both actual and predicted relative reductions in the sum of squares
3058 ! are at most FTOL.
3059 ! 2, relative error between two consecutive iterates is at most XTOL.
3060 ! 3, conditions for INFO = 1 and INFO = 2 both hold.
3061 ! 4, the cosine of the angle between FVEC and any column of the jacobian
3062 ! is at most GTOL in absolute value.
3063 ! 5, number of calls to FCN has reached or exceeded MAXFEV.
3064 ! 6, FTOL is too small. No further reduction in the sum of squares
3065 ! is possible.
3066 ! 7, XTOL is too small. No further improvement in the approximate
3067 ! solution X is possible.
3068 ! 8, GTOL is too small. FVEC is orthogonal to the columns of the
3069 ! jacobian to machine precision.
3070 !
3071 ! Output, integer ( kind = 4 ) NFEV, the number of calls to FCN.
3072 !
3073 ! Output, real ( kind = 8 ) FJAC(LDFJAC,N), an M by N array. The upper
3074 ! N by N submatrix of FJAC contains an upper triangular matrix R with
3075 ! diagonal elements of nonincreasing magnitude such that
3076 !
3077 ! P' * ( JAC' * JAC ) * P = R' * R,
3078 !
3079 ! where P is a permutation matrix and JAC is the final calculated jacobian.
3080 ! Column J of P is column IPVT(J) of the identity matrix. The lower
3081 ! trapezoidal part of FJAC contains information generated during
3082 ! the computation of R.
3083 !
3084 ! Input, integer ( kind = 4 ) LDFJAC, the leading dimension of FJAC.
3085 ! LDFJAC must be at least M.
3086 !
3087 ! Output, integer ( kind = 4 ) IPVT(N), defines a permutation matrix P such
3088 ! that JAC * P = Q * R, where JAC is the final calculated jacobian, Q is
3089 ! orthogonal (not stored), and R is upper triangular with diagonal
3090 ! elements of nonincreasing magnitude. Column J of P is column IPVT(J)
3091 ! of the identity matrix.
3092 !
3093 ! Output, real ( kind = 8 ) QTF(N), the first N elements of Q'*FVEC.
3094 !
3095  implicit none
3096 
3097  integer(kind=4) ldfjac
3098  integer(kind=4) m
3099  integer(kind=4) n
3100 
3101  real(kind=8) actred
3102  real(kind=8) delta
3103  real(kind=8) diag(n)
3104  real(kind=8) dirder
3105  real(kind=8) enorm
3106  real(kind=8) epsfcn
3107  real(kind=8) epsmch
3108  real(kind=8) factor
3109  external fcn
3110  real(kind=8) fjac(ldfjac, n)
3111  real(kind=8) fnorm
3112  real(kind=8) fnorm1
3113  real(kind=8) ftol
3114  real(kind=8) :: fvec(m), xdat(m), ydat(m)
3115  real(kind=8) gnorm
3116  real(kind=8) gtol
3117  integer(kind=4) i
3118  integer(kind=4) iflag
3119  integer(kind=4) iter
3120  integer(kind=4) info
3121  integer(kind=4) ipvt(n)
3122  integer(kind=4) j
3123  integer(kind=4) l
3124  integer(kind=4) maxfev
3125  integer(kind=4) mode
3126  integer(kind=4) nfev
3127  integer(kind=4) nprint
3128  real(kind=8) par
3129  logical pivot
3130  real(kind=8) pnorm
3131  real(kind=8) prered
3132  real(kind=8) qtf(n)
3133  real(kind=8) ratio
3134  real(kind=8) sum2
3135  real(kind=8) temp
3136  real(kind=8) temp1
3137  real(kind=8) temp2
3138  real(kind=8) wa1(n)
3139  real(kind=8) wa2(n)
3140  real(kind=8) wa3(n)
3141  real(kind=8) wa4(m)
3142  real(kind=8) x(n)
3143  real(kind=8) xnorm
3144  real(kind=8) xtol
3145 
3146  epsmch = epsilon(epsmch)
3147 
3148  info = 0
3149  iflag = 0
3150  nfev = 0
3151 
3152  if (n <= 0) then
3153  go to 300
3154  else if (m < n) then
3155  go to 300
3156  else if (ldfjac < m) then
3157  go to 300
3158  else if (ftol < 0.0d+00) then
3159  go to 300
3160  else if (xtol < 0.0d+00) then
3161  go to 300
3162  else if (gtol < 0.0d+00) then
3163  go to 300
3164  else if (maxfev <= 0) then
3165  go to 300
3166  else if (factor <= 0.0d+00) then
3167  go to 300
3168  end if
3169 
3170  if (mode == 2) then
3171  do j = 1, n
3172  if (diag(j) <= 0.0d+00) then
3173  go to 300
3174  end if
3175  end do
3176  end if
3177 !
3178 ! Evaluate the function at the starting point and calculate its norm.
3179 !
3180  iflag = 1
3181  call fcn(m, n, x, xdat, ydat, fvec, iflag)
3182  nfev = 1
3183 
3184  if (iflag < 0) then
3185  go to 300
3186  end if
3187 
3188  fnorm = enorm(m, fvec)
3189 !
3190 ! Initialize Levenberg-Marquardt parameter and iteration counter.
3191 !
3192  par = 0.0d+00
3193  iter = 1
3194 !
3195 ! Beginning of the outer loop.
3196 !
3197 30 continue
3198 !
3199 ! Calculate the jacobian matrix.
3200 !
3201  iflag = 2
3202  call fdjac2(fcn, m, n, x, xdat, ydat, fvec, fjac, ldfjac, iflag, epsfcn)
3203  nfev = nfev + n
3204 
3205  if (iflag < 0) then
3206  go to 300
3207  end if
3208 !
3209 ! If requested, call FCN to enable printing of iterates.
3210 !
3211  if (0 < nprint) then
3212  iflag = 0
3213  if (mod(iter - 1, nprint) == 0) then
3214  call fcn(m, n, x, xdat, ydat, fvec, iflag)
3215  end if
3216  if (iflag < 0) then
3217  go to 300
3218  end if
3219  end if
3220 !
3221 ! Compute the QR factorization of the jacobian.
3222 !
3223  pivot = .true.
3224  call qrfac(m, n, fjac, ldfjac, pivot, ipvt, n, wa1, wa2)
3225 !
3226 ! On the first iteration and if MODE is 1, scale according
3227 ! to the norms of the columns of the initial jacobian.
3228 !
3229  if (iter == 1) then
3230 
3231  if (mode /= 2) then
3232  diag(1:n) = wa2(1:n)
3233  do j = 1, n
3234  if (wa2(j) == 0.0d+00) then
3235  diag(j) = 1.0d+00
3236  end if
3237  end do
3238  end if
3239 !
3240 ! On the first iteration, calculate the norm of the scaled X
3241 ! and initialize the step bound DELTA.
3242 !
3243  wa3(1:n) = diag(1:n)*x(1:n)
3244  xnorm = enorm(n, wa3)
3245  delta = factor*xnorm
3246  if (delta == 0.0d+00) then
3247  delta = factor
3248  end if
3249  end if
3250 !
3251 ! Form Q' * FVEC and store the first N components in QTF.
3252 !
3253  wa4(1:m) = fvec(1:m)
3254 
3255  do j = 1, n
3256 
3257  if (fjac(j, j) /= 0.0d+00) then
3258  sum2 = dot_product(wa4(j:m), fjac(j:m, j))
3259  temp = -sum2/fjac(j, j)
3260  wa4(j:m) = wa4(j:m) + fjac(j:m, j)*temp
3261  end if
3262 
3263  fjac(j, j) = wa1(j)
3264  qtf(j) = wa4(j)
3265 
3266  end do
3267 !
3268 ! Compute the norm of the scaled gradient.
3269 !
3270  gnorm = 0.0d+00
3271 
3272  if (fnorm /= 0.0d+00) then
3273 
3274  do j = 1, n
3275 
3276  l = ipvt(j)
3277 
3278  if (wa2(l) /= 0.0d+00) then
3279  sum2 = 0.0d+00
3280  do i = 1, j
3281  sum2 = sum2 + fjac(i, j)*(qtf(i)/fnorm)
3282  end do
3283  gnorm = max(gnorm, abs(sum2/wa2(l)))
3284  end if
3285 
3286  end do
3287 
3288  end if
3289 !
3290 ! Test for convergence of the gradient norm.
3291 !
3292  if (gnorm <= gtol) then
3293  info = 4
3294  go to 300
3295  end if
3296 !
3297 ! Rescale if necessary.
3298 !
3299  if (mode /= 2) then
3300  do j = 1, n
3301  diag(j) = max(diag(j), wa2(j))
3302  end do
3303  end if
3304 !
3305 ! Beginning of the inner loop.
3306 !
3307 200 continue
3308 !
3309 ! Determine the Levenberg-Marquardt parameter.
3310 !
3311  call lmpar(n, fjac, ldfjac, ipvt, diag, qtf, delta, par, wa1, wa2)
3312 !
3313 ! Store the direction P and X + P.
3314 ! Calculate the norm of P.
3315 !
3316  wa1(1:n) = -wa1(1:n)
3317  wa2(1:n) = x(1:n) + wa1(1:n)
3318  wa3(1:n) = diag(1:n)*wa1(1:n)
3319 
3320  pnorm = enorm(n, wa3)
3321 !
3322 ! On the first iteration, adjust the initial step bound.
3323 !
3324  if (iter == 1) then
3325  delta = min(delta, pnorm)
3326  end if
3327 !
3328 ! Evaluate the function at X + P and calculate its norm.
3329 !
3330  iflag = 1
3331  call fcn(m, n, wa2, xdat, ydat, wa4, iflag)
3332  nfev = nfev + 1
3333  if (iflag < 0) then
3334  go to 300
3335  end if
3336  fnorm1 = enorm(m, wa4)
3337 !
3338 ! Compute the scaled actual reduction.
3339 !
3340  if (0.1d+00*fnorm1 < fnorm) then
3341  actred = 1.0d+00 - (fnorm1/fnorm)**2
3342  else
3343  actred = -1.0d+00
3344  end if
3345 !
3346 ! Compute the scaled predicted reduction and the scaled directional derivative.
3347 !
3348  do j = 1, n
3349  wa3(j) = 0.0d+00
3350  l = ipvt(j)
3351  temp = wa1(l)
3352  wa3(1:j) = wa3(1:j) + fjac(1:j, j)*temp
3353  end do
3354 
3355  temp1 = enorm(n, wa3)/fnorm
3356  temp2 = (sqrt(par)*pnorm)/fnorm
3357  prered = temp1**2 + temp2**2/0.5d+00
3358  dirder = -(temp1**2 + temp2**2)
3359 !
3360 ! Compute the ratio of the actual to the predicted reduction.
3361 !
3362  ratio = 0.0d+00
3363  if (prered /= 0.0d+00) then
3364  ratio = actred/prered
3365  end if
3366 !
3367 ! Update the step bound.
3368 !
3369  if (ratio <= 0.25d+00) then
3370 
3371  if (actred >= 0.0d+00) then
3372  temp = 0.5d+00
3373  endif
3374 
3375  if (actred < 0.0d+00) then
3376  temp = 0.5d+00*dirder/(dirder + 0.5d+00*actred)
3377  end if
3378 
3379  if (0.1d+00*fnorm1 >= fnorm .or. temp < 0.1d+00) then
3380  temp = 0.1d+00
3381  end if
3382 
3383  delta = temp*min(delta, pnorm/0.1d+00)
3384  par = par/temp
3385 
3386  else
3387 
3388  if (par == 0.0d+00 .or. ratio >= 0.75d+00) then
3389  delta = 2.0d+00*pnorm
3390  par = 0.5d+00*par
3391  end if
3392 
3393  end if
3394 !
3395 ! Test for successful iteration.
3396 !
3397 
3398 !
3399 ! Successful iteration. update X, FVEC, and their norms.
3400 !
3401  if (0.0001d+00 <= ratio) then
3402  x(1:n) = wa2(1:n)
3403  wa2(1:n) = diag(1:n)*x(1:n)
3404  fvec(1:m) = wa4(1:m)
3405  xnorm = enorm(n, wa2)
3406  fnorm = fnorm1
3407  iter = iter + 1
3408  end if
3409 !
3410 ! Tests for convergence.
3411 !
3412  if (abs(actred) <= ftol .and. prered <= ftol &
3413  .and. 0.5d+00*ratio <= 1.0d+00) then
3414  info = 1
3415  end if
3416 
3417  if (delta <= xtol*xnorm) then
3418  info = 2
3419  end if
3420 
3421  if (abs(actred) <= ftol .and. prered <= ftol &
3422  .and. 0.5d+00*ratio <= 1.0d+00 .and. info == 2) info = 3
3423 
3424  if (info /= 0) then
3425  go to 300
3426  end if
3427 !
3428 ! Tests for termination and stringent tolerances.
3429 !
3430  if (maxfev <= nfev) then
3431  info = 5
3432  end if
3433 
3434  if (abs(actred) <= epsmch .and. prered <= epsmch &
3435  .and. 0.5d+00*ratio <= 1.0d+00) then
3436  info = 6
3437  end if
3438 
3439  if (delta <= epsmch*xnorm) then
3440  info = 7
3441  end if
3442 
3443  if (gnorm <= epsmch) then
3444  info = 8
3445  end if
3446 
3447  if (info /= 0) then
3448  go to 300
3449  end if
3450 !
3451 ! End of the inner loop. Repeat if iteration unsuccessful.
3452 !
3453  if (ratio < 0.0001d+00) then
3454  go to 200
3455  end if
3456 !
3457 ! End of the outer loop.
3458 !
3459  go to 30
3460 
3461 300 continue
3462 !
3463 ! Termination, either normal or user imposed.
3464 !
3465  if (iflag < 0) then
3466  info = iflag
3467  end if
3468 
3469  iflag = 0
3470 
3471  if (0 < nprint) then
3472  call fcn(m, n, x, xdat, ydat, fvec, iflag)
3473  end if
3474 
3475  return
real(kind=8) function enorm(n, x)
subroutine qrfac(m, n, a, lda, pivot, ipvt, lipvt, rdiag, acnorm)
subroutine lmpar(n, r, ldr, ipvt, diag, qtb, delta, par, x, sdiag)
subroutine fdjac2(fcn, m, n, x, xdat, ydat, fvec, fjac, ldfjac, iflag, epsfcn)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ lmdif1()

subroutine lmdif1 ( external  fcn,
integer(kind=4)  m,
integer(kind=4)  n,
real(kind=8), dimension(n)  x,
real(kind=8), dimension(m)  xdat,
real(kind=8), dimension(m)  ydat,
real(kind=8), dimension(m)  fvec,
real(kind=8)  tol,
integer(kind=4)  info 
)

Definition at line 3478 of file suews_util_minpack.f95.

References lmdif().

Referenced by anohm_module::anohm_shapefit().

3478 
3479 !*****************************************************************************80
3480 !
3481 !! LMDIF1 minimizes M functions in N variables using Levenberg-Marquardt method.
3482 !
3483 ! Discussion:
3484 !
3485 ! LMDIF1 minimizes the sum of the squares of M nonlinear functions in
3486 ! N variables by a modification of the Levenberg-Marquardt algorithm.
3487 ! This is done by using the more general least-squares solver LMDIF.
3488 ! The user must provide a subroutine which calculates the functions.
3489 ! The jacobian is then calculated by a forward-difference approximation.
3490 !
3491 ! Licensing:
3492 !
3493 ! This code is distributed under the GNU LGPL license.
3494 !
3495 ! Modified:
3496 !
3497 ! 20 Jul 2017
3498 !
3499 ! Author:
3500 !
3501 ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom.
3502 ! FORTRAN90 version by John Burkardt.
3503 ! adapted for AnOHM in SUEWS, Ting Sun, ting.sun@reading.ac.uk
3504 !
3505 ! Reference:
3506 !
3507 ! Jorge More, Burton Garbow, Kenneth Hillstrom,
3508 ! User Guide for MINPACK-1,
3509 ! Technical Report ANL-80-74,
3510 ! Argonne National Laboratory, 1980.
3511 !
3512 ! Parameters:
3513 !
3514 ! Input, external FCN, the name of the user-supplied subroutine which
3515 ! calculates the functions. The routine should have the form:
3516 ! subroutine fcn ( m, n, x, xdat, ydat, fvec, iflag ) ! xdat, ydat added for AnOHM, TS 20 Jul 2017
3517 ! integer ( kind = 4 ) m
3518 ! integer ( kind = 4 ) n
3519 ! real ( kind = 8 ) fvec(m),xdat(m),ydat(m)
3520 ! integer ( kind = 4 ) iflag
3521 ! real ( kind = 8 ) x(n)
3522 ! If IFLAG = 0 on input, then FCN is only being called to allow the user
3523 ! to print out the current iterate.
3524 ! If IFLAG = 1 on input, FCN should calculate the functions at X and
3525 ! return this vector in FVEC.
3526 ! To terminate the algorithm, FCN may set IFLAG negative on return.
3527 !
3528 ! Input, integer ( kind = 4 ) M, the number of functions.
3529 !
3530 ! Input, integer ( kind = 4 ) N, the number of variables.
3531 ! N must not exceed M.
3532 !
3533 ! Input/output, real ( kind = 8 ) X(N). On input, X must contain an initial
3534 ! estimate of the solution vector. On output X contains the final
3535 ! estimate of the solution vector.
3536 
3537 ! Input, real ( kind = 8 ) XDAT(M), YDAT(M). On input, XDAT/YDAT must contain
3538 ! observations to construct FVEC(M). ! added for AnOHM, TS 20 Jul 2017
3539 !
3540 ! Output, real ( kind = 8 ) FVEC(M), the functions evaluated at the output X.
3541 !
3542 ! Input, real ( kind = 8 ) TOL. Termination occurs when the algorithm
3543 ! estimates either that the relative error in the sum of squares is at
3544 ! most TOL or that the relative error between X and the solution is at
3545 ! most TOL. TOL should be nonnegative.
3546 !
3547 ! Output, integer ( kind = 4 ) INFO, error flag. If the user has terminated
3548 ! execution, INFO is set to the (negative) value of IFLAG. See description
3549 ! of FCN. Otherwise, INFO is set as follows:
3550 ! 0, improper input parameters.
3551 ! 1, algorithm estimates that the relative error in the sum of squares
3552 ! is at most TOL.
3553 ! 2, algorithm estimates that the relative error between X and the
3554 ! solution is at most TOL.
3555 ! 3, conditions for INFO = 1 and INFO = 2 both hold.
3556 ! 4, FVEC is orthogonal to the columns of the jacobian to machine precision.
3557 ! 5, number of calls to FCN has reached or exceeded 200*(N+1).
3558 ! 6, TOL is too small. No further reduction in the sum of squares
3559 ! is possible.
3560 ! 7, TOL is too small. No further improvement in the approximate
3561 ! solution X is possible.
3562 !
3563  implicit none
3564 
3565  integer(kind=4) m
3566  integer(kind=4) n
3567 
3568  real(kind=8) diag(n)
3569  real(kind=8) epsfcn
3570  real(kind=8) factor
3571  external fcn
3572  real(kind=8) fjac(m, n)
3573  real(kind=8) ftol
3574  real(kind=8) ::fvec(m), xdat(m), ydat(m)
3575  real(kind=8) gtol
3576  integer(kind=4) info
3577  integer(kind=4) ipvt(n)
3578  integer(kind=4) ldfjac
3579  integer(kind=4) maxfev
3580  integer(kind=4) mode
3581  integer(kind=4) nfev
3582  integer(kind=4) nprint
3583  real(kind=8) qtf(n)
3584  real(kind=8) tol
3585  real(kind=8) x(n)
3586  real(kind=8) xtol
3587 
3588  info = 0
3589 
3590  if (n <= 0) then
3591  return
3592  else if (m < n) then
3593  return
3594  else if (tol < 0.0d+00) then
3595  return
3596  end if
3597 
3598  factor = 100.0d+00
3599  maxfev = 200*(n + 1)
3600  ftol = tol
3601  xtol = tol
3602  gtol = 0.0d+00
3603  epsfcn = 0.0d+00
3604  mode = 1
3605  nprint = 0
3606  ldfjac = m
3607 ! print*, 'x in limdif1',x
3608  call lmdif(fcn, m, n, x, xdat, ydat, fvec, ftol, xtol, gtol, maxfev, epsfcn, &
3609  diag, mode, factor, nprint, info, nfev, fjac, ldfjac, ipvt, qtf)
3610 
3611  if (info == 8) then
3612  info = 4
3613  end if
3614 
3615  return
subroutine lmdif(fcn, m, n, x, xdat, ydat, fvec, ftol, xtol, gtol, maxfev, epsfcn, diag, mode, factor, nprint, info, nfev, fjac, ldfjac, ipvt, qtf)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ lmpar()

subroutine lmpar ( integer(kind=4)  n,
real(kind=8), dimension(ldr, n)  r,
integer(kind=4)  ldr,
integer(kind=4), dimension(n)  ipvt,
real(kind=8), dimension(n)  diag,
real(kind=8), dimension(n)  qtb,
real(kind=8)  delta,
real(kind=8)  par,
real(kind=8), dimension(n)  x,
real(kind=8), dimension(n)  sdiag 
)

Definition at line 3618 of file suews_util_minpack.f95.

References qrsolv().

Referenced by lmder(), lmdif(), and lmstr().

3618 
3619 !*****************************************************************************80
3620 !
3621 !! LMPAR computes a parameter for the Levenberg-Marquardt method.
3622 !
3623 ! Discussion:
3624 !
3625 ! Given an M by N matrix A, an N by N nonsingular diagonal
3626 ! matrix D, an M-vector B, and a positive number DELTA,
3627 ! the problem is to determine a value for the parameter
3628 ! PAR such that if X solves the system
3629 !
3630 ! A*X = B,
3631 ! sqrt ( PAR ) * D * X = 0,
3632 !
3633 ! in the least squares sense, and DXNORM is the euclidean
3634 ! norm of D*X, then either PAR is zero and
3635 !
3636 ! ( DXNORM - DELTA ) <= 0.1 * DELTA,
3637 !
3638 ! or PAR is positive and
3639 !
3640 ! abs ( DXNORM - DELTA) <= 0.1 * DELTA.
3641 !
3642 ! This subroutine completes the solution of the problem
3643 ! if it is provided with the necessary information from the
3644 ! QR factorization, with column pivoting, of A. That is, if
3645 ! A*P = Q*R, where P is a permutation matrix, Q has orthogonal
3646 ! columns, and R is an upper triangular matrix with diagonal
3647 ! elements of nonincreasing magnitude, then LMPAR expects
3648 ! the full upper triangle of R, the permutation matrix P,
3649 ! and the first N components of Q'*B. On output
3650 ! LMPAR also provides an upper triangular matrix S such that
3651 !
3652 ! P' * ( A' * A + PAR * D * D ) * P = S'* S.
3653 !
3654 ! S is employed within LMPAR and may be of separate interest.
3655 !
3656 ! Only a few iterations are generally needed for convergence
3657 ! of the algorithm. If, however, the limit of 10 iterations
3658 ! is reached, then the output PAR will contain the best
3659 ! value obtained so far.
3660 !
3661 ! Licensing:
3662 !
3663 ! This code is distributed under the GNU LGPL license.
3664 !
3665 ! Modified:
3666 !
3667 ! 24 January 2014
3668 !
3669 ! Author:
3670 !
3671 ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom.
3672 ! FORTRAN90 version by John Burkardt.
3673 !
3674 ! Reference:
3675 !
3676 ! Jorge More, Burton Garbow, Kenneth Hillstrom,
3677 ! User Guide for MINPACK-1,
3678 ! Technical Report ANL-80-74,
3679 ! Argonne National Laboratory, 1980.
3680 !
3681 ! Parameters:
3682 !
3683 ! Input, integer ( kind = 4 ) N, the order of R.
3684 !
3685 ! Input/output, real ( kind = 8 ) R(LDR,N),the N by N matrix. The full
3686 ! upper triangle must contain the full upper triangle of the matrix R.
3687 ! On output the full upper triangle is unaltered, and the strict lower
3688 ! triangle contains the strict upper triangle (transposed) of the upper
3689 ! triangular matrix S.
3690 !
3691 ! Input, integer ( kind = 4 ) LDR, the leading dimension of R. LDR must be
3692 ! no less than N.
3693 !
3694 ! Input, integer ( kind = 4 ) IPVT(N), defines the permutation matrix P
3695 ! such that A*P = Q*R. Column J of P is column IPVT(J) of the
3696 ! identity matrix.
3697 !
3698 ! Input, real ( kind = 8 ) DIAG(N), the diagonal elements of the matrix D.
3699 !
3700 ! Input, real ( kind = 8 ) QTB(N), the first N elements of the vector Q'*B.
3701 !
3702 ! Input, real ( kind = 8 ) DELTA, an upper bound on the euclidean norm
3703 ! of D*X. DELTA should be positive.
3704 !
3705 ! Input/output, real ( kind = 8 ) PAR. On input an initial estimate of the
3706 ! Levenberg-Marquardt parameter. On output the final estimate.
3707 ! PAR should be nonnegative.
3708 !
3709 ! Output, real ( kind = 8 ) X(N), the least squares solution of the system
3710 ! A*X = B, sqrt(PAR)*D*X = 0, for the output value of PAR.
3711 !
3712 ! Output, real ( kind = 8 ) SDIAG(N), the diagonal elements of the upper
3713 ! triangular matrix S.
3714 !
3715  implicit none
3716 
3717  integer(kind=4) ldr
3718  integer(kind=4) n
3719 
3720  real(kind=8) delta
3721  real(kind=8) diag(n)
3722  real(kind=8) dwarf
3723  real(kind=8) dxnorm
3724  real(kind=8) enorm
3725  real(kind=8) gnorm
3726  real(kind=8) fp
3727  ! integer ( kind = 4 ) i
3728  integer(kind=4) ipvt(n)
3729  integer(kind=4) iter
3730  integer(kind=4) j
3731  integer(kind=4) k
3732  integer(kind=4) l
3733  integer(kind=4) nsing
3734  real(kind=8) par
3735  real(kind=8) parc
3736  real(kind=8) parl
3737  real(kind=8) paru
3738  ! real ( kind = 8 ) qnorm
3739  real(kind=8) qtb(n)
3740  real(kind=8) r(ldr, n)
3741  real(kind=8) sdiag(n)
3742  real(kind=8) sum2
3743  real(kind=8) temp
3744  real(kind=8) wa1(n)
3745  real(kind=8) wa2(n)
3746  real(kind=8) x(n)
3747 !
3748 ! DWARF is the smallest positive magnitude.
3749 !
3750  dwarf = tiny(dwarf)
3751 !
3752 ! Compute and store in X the Gauss-Newton direction.
3753 !
3754 ! If the jacobian is rank-deficient, obtain a least squares solution.
3755 !
3756  nsing = n
3757 
3758  do j = 1, n
3759  wa1(j) = qtb(j)
3760  if (r(j, j) == 0.0d+00 .and. nsing == n) then
3761  nsing = j - 1
3762  end if
3763  if (nsing < n) then
3764  wa1(j) = 0.0d+00
3765  end if
3766  end do
3767 
3768  do k = 1, nsing
3769  j = nsing - k + 1
3770  wa1(j) = wa1(j)/r(j, j)
3771  temp = wa1(j)
3772  wa1(1:j - 1) = wa1(1:j - 1) - r(1:j - 1, j)*temp
3773  end do
3774 
3775  do j = 1, n
3776  l = ipvt(j)
3777  x(l) = wa1(j)
3778  end do
3779 !
3780 ! Initialize the iteration counter.
3781 ! Evaluate the function at the origin, and test
3782 ! for acceptance of the Gauss-Newton direction.
3783 !
3784  iter = 0
3785  wa2(1:n) = diag(1:n)*x(1:n)
3786  dxnorm = enorm(n, wa2)
3787  fp = dxnorm - delta
3788 
3789  if (fp <= 0.1d+00*delta) then
3790  if (iter == 0) then
3791  par = 0.0d+00
3792  end if
3793  return
3794  end if
3795 !
3796 ! If the jacobian is not rank deficient, the Newton
3797 ! step provides a lower bound, PARL, for the zero of
3798 ! the function.
3799 !
3800 ! Otherwise set this bound to zero.
3801 !
3802  parl = 0.0d+00
3803 
3804  if (n <= nsing) then
3805 
3806  do j = 1, n
3807  l = ipvt(j)
3808  wa1(j) = diag(l)*(wa2(l)/dxnorm)
3809  end do
3810 
3811  do j = 1, n
3812  sum2 = dot_product(wa1(1:j - 1), r(1:j - 1, j))
3813  wa1(j) = (wa1(j) - sum2)/r(j, j)
3814  end do
3815 
3816  temp = enorm(n, wa1)
3817  parl = ((fp/delta)/temp)/temp
3818 
3819  end if
3820 !
3821 ! Calculate an upper bound, PARU, for the zero of the function.
3822 !
3823  do j = 1, n
3824  sum2 = dot_product(qtb(1:j), r(1:j, j))
3825  l = ipvt(j)
3826  wa1(j) = sum2/diag(l)
3827  end do
3828 
3829  gnorm = enorm(n, wa1)
3830  paru = gnorm/delta
3831 
3832  if (paru == 0.0d+00) then
3833  paru = dwarf/min(delta, 0.1d+00)
3834  end if
3835 !
3836 ! If the input PAR lies outside of the interval (PARL, PARU),
3837 ! set PAR to the closer endpoint.
3838 !
3839  par = max(par, parl)
3840  par = min(par, paru)
3841  if (par == 0.0d+00) then
3842  par = gnorm/dxnorm
3843  end if
3844 !
3845 ! Beginning of an iteration.
3846 !
3847  do
3848 
3849  iter = iter + 1
3850 !
3851 ! Evaluate the function at the current value of PAR.
3852 !
3853  if (par == 0.0d+00) then
3854  par = max(dwarf, 0.001d+00*paru)
3855  end if
3856 
3857  wa1(1:n) = sqrt(par)*diag(1:n)
3858 
3859  call qrsolv(n, r, ldr, ipvt, wa1, qtb, x, sdiag)
3860 
3861  wa2(1:n) = diag(1:n)*x(1:n)
3862  dxnorm = enorm(n, wa2)
3863  temp = fp
3864  fp = dxnorm - delta
3865 !
3866 ! If the function is small enough, accept the current value of PAR.
3867 !
3868  if (abs(fp) <= 0.1d+00*delta) then
3869  exit
3870  end if
3871 !
3872 ! Test for the exceptional cases where PARL
3873 ! is zero or the number of iterations has reached 10.
3874 !
3875  if (parl == 0.0d+00 .and. fp <= temp .and. temp < 0.0d+00) then
3876  exit
3877  else if (iter == 10) then
3878  exit
3879  end if
3880 !
3881 ! Compute the Newton correction.
3882 !
3883  do j = 1, n
3884  l = ipvt(j)
3885  wa1(j) = diag(l)*(wa2(l)/dxnorm)
3886  end do
3887 
3888  do j = 1, n
3889  wa1(j) = wa1(j)/sdiag(j)
3890  temp = wa1(j)
3891  wa1(j + 1:n) = wa1(j + 1:n) - r(j + 1:n, j)*temp
3892  end do
3893 
3894  temp = enorm(n, wa1)
3895  parc = ((fp/delta)/temp)/temp
3896 !
3897 ! Depending on the sign of the function, update PARL or PARU.
3898 !
3899  if (0.0d+00 < fp) then
3900  parl = max(parl, par)
3901  else if (fp < 0.0d+00) then
3902  paru = min(paru, par)
3903  end if
3904 !
3905 ! Compute an improved estimate for PAR.
3906 !
3907  par = max(parl, par + parc)
3908 !
3909 ! End of an iteration.
3910 !
3911  end do
3912 !
3913 ! Termination.
3914 !
3915  if (iter == 0) then
3916  par = 0.0d+00
3917  end if
3918 
3919  return
subroutine qrsolv(n, r, ldr, ipvt, diag, qtb, x, sdiag)
real(kind=8) function enorm(n, x)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ lmstr()

subroutine lmstr ( external  fcn,
integer(kind=4)  m,
integer(kind=4)  n,
real(kind=8), dimension(n)  x,
real(kind=8), dimension(m)  fvec,
real(kind=8), dimension(ldfjac, n)  fjac,
integer(kind=4)  ldfjac,
real(kind=8)  ftol,
real(kind=8)  xtol,
real(kind=8)  gtol,
integer(kind=4)  maxfev,
real(kind=8), dimension(n)  diag,
integer(kind=4)  mode,
real(kind=8)  factor,
integer(kind=4)  nprint,
integer(kind=4)  info,
integer(kind=4)  nfev,
integer(kind=4)  njev,
integer(kind=4), dimension(n)  ipvt,
real(kind=8), dimension(n)  qtf 
)

Definition at line 3923 of file suews_util_minpack.f95.

References lmpar(), qrfac(), and rwupdt().

Referenced by lmstr1().

3923 
3924 !*****************************************************************************80
3925 !
3926 !! LMSTR minimizes M functions in N variables using Levenberg-Marquardt method.
3927 !
3928 ! Discussion:
3929 !
3930 ! LMSTR minimizes the sum of the squares of M nonlinear functions in
3931 ! N variables by a modification of the Levenberg-Marquardt algorithm
3932 ! which uses minimal storage.
3933 !
3934 ! The user must provide a subroutine which calculates the functions and
3935 ! the rows of the jacobian.
3936 !
3937 ! Licensing:
3938 !
3939 ! This code is distributed under the GNU LGPL license.
3940 !
3941 ! Modified:
3942 !
3943 ! 06 April 2010
3944 !
3945 ! Author:
3946 !
3947 ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom.
3948 ! FORTRAN90 version by John Burkardt.
3949 !
3950 ! Reference:
3951 !
3952 ! Jorge More, Burton Garbow, Kenneth Hillstrom,
3953 ! User Guide for MINPACK-1,
3954 ! Technical Report ANL-80-74,
3955 ! Argonne National Laboratory, 1980.
3956 !
3957 ! Parameters:
3958 !
3959 ! Input, external FCN, the name of the user-supplied subroutine which
3960 ! calculates the functions and the rows of the jacobian.
3961 ! FCN should have the form:
3962 ! subroutine fcn ( m, n, x, fvec, fjrow, iflag )
3963 ! integer ( kind = 4 ) m
3964 ! integer ( kind = 4 ) n
3965 ! real ( kind = 8 ) fjrow(n)
3966 ! real ( kind = 8 ) fvec(m)
3967 ! integer ( kind = 4 ) iflag
3968 ! real ( kind = 8 ) x(n)
3969 ! If IFLAG = 0 on input, then FCN is only being called to allow the user
3970 ! to print out the current iterate.
3971 ! If IFLAG = 1 on input, FCN should calculate the functions at X and
3972 ! return this vector in FVEC.
3973 ! If the input value of IFLAG is I > 1, calculate the (I-1)-st row of
3974 ! the jacobian at X, and return this vector in FJROW.
3975 ! To terminate the algorithm, set the output value of IFLAG negative.
3976 !
3977 ! Input, integer ( kind = 4 ) M, the number of functions.
3978 !
3979 ! Input, integer ( kind = 4 ) N, the number of variables.
3980 ! N must not exceed M.
3981 !
3982 ! Input/output, real ( kind = 8 ) X(N). On input, X must contain an initial
3983 ! estimate of the solution vector. On output X contains the final
3984 ! estimate of the solution vector.
3985 !
3986 ! Output, real ( kind = 8 ) FVEC(M), the functions evaluated at the output X.
3987 !
3988 ! Output, real ( kind = 8 ) FJAC(LDFJAC,N), an N by N array. The upper
3989 ! triangle of FJAC contains an upper triangular matrix R such that
3990 !
3991 ! P' * ( JAC' * JAC ) * P = R' * R,
3992 !
3993 ! where P is a permutation matrix and JAC is the final calculated jacobian.
3994 ! Column J of P is column IPVT(J) of the identity matrix. The lower
3995 ! triangular part of FJAC contains information generated during
3996 ! the computation of R.
3997 !
3998 ! Input, integer ( kind = 4 ) LDFJAC, the leading dimension of FJAC.
3999 ! LDFJAC must be at least N.
4000 !
4001 ! Input, real ( kind = 8 ) FTOL. Termination occurs when both the actual and
4002 ! predicted relative reductions in the sum of squares are at most FTOL.
4003 ! Therefore, FTOL measures the relative error desired in the sum of
4004 ! squares. FTOL should be nonnegative.
4005 !
4006 ! Input, real ( kind = 8 ) XTOL. Termination occurs when the relative error
4007 ! between two consecutive iterates is at most XTOL. XTOL should be
4008 ! nonnegative.
4009 !
4010 ! Input, real ( kind = 8 ) GTOL. termination occurs when the cosine of the
4011 ! angle between FVEC and any column of the jacobian is at most GTOL in
4012 ! absolute value. Therefore, GTOL measures the orthogonality desired
4013 ! between the function vector and the columns of the jacobian. GTOL should
4014 ! be nonnegative.
4015 !
4016 ! Input, integer ( kind = 4 ) MAXFEV. Termination occurs when the number
4017 ! of calls to FCN with IFLAG = 1 is at least MAXFEV by the end of
4018 ! an iteration.
4019 !
4020 ! Input/output, real ( kind = 8 ) DIAG(N). If MODE = 1, then DIAG is set
4021 ! internally. If MODE = 2, then DIAG must contain positive entries that
4022 ! serve as multiplicative scale factors for the variables.
4023 !
4024 ! Input, integer ( kind = 4 ) MODE, scaling option.
4025 ! 1, variables will be scaled internally.
4026 ! 2, scaling is specified by the input DIAG vector.
4027 !
4028 ! Input, real ( kind = 8 ) FACTOR, determines the initial step bound. This
4029 ! bound is set to the product of FACTOR and the euclidean norm of DIAG*X if
4030 ! nonzero, or else to FACTOR itself. In most cases, FACTOR should lie
4031 ! in the interval (0.1, 100) with 100 the recommended value.
4032 !
4033 ! Input, integer ( kind = 4 ) NPRINT, enables controlled printing of iterates
4034 ! if it is positive. In this case, FCN is called with IFLAG = 0 at the
4035 ! beginning of the first iteration and every NPRINT iterations thereafter
4036 ! and immediately prior to return, with X and FVEC available
4037 ! for printing. If NPRINT is not positive, no special calls
4038 ! of FCN with IFLAG = 0 are made.
4039 !
4040 ! Output, integer ( kind = 4 ) INFO, error flag. If the user has terminated
4041 ! execution, INFO is set to the (negative) value of IFLAG. See the
4042 ! description of FCN. Otherwise, INFO is set as follows:
4043 ! 0, improper input parameters.
4044 ! 1, both actual and predicted relative reductions in the sum of squares
4045 ! are at most FTOL.
4046 ! 2, relative error between two consecutive iterates is at most XTOL.
4047 ! 3, conditions for INFO = 1 and INFO = 2 both hold.
4048 ! 4, the cosine of the angle between FVEC and any column of the jacobian
4049 ! is at most GTOL in absolute value.
4050 ! 5, number of calls to FCN with IFLAG = 1 has reached MAXFEV.
4051 ! 6, FTOL is too small. No further reduction in the sum of squares is
4052 ! possible.
4053 ! 7, XTOL is too small. No further improvement in the approximate
4054 ! solution X is possible.
4055 ! 8, GTOL is too small. FVEC is orthogonal to the columns of the
4056 ! jacobian to machine precision.
4057 !
4058 ! Output, integer ( kind = 4 ) NFEV, the number of calls to FCN
4059 ! with IFLAG = 1.
4060 !
4061 ! Output, integer ( kind = 4 ) NJEV, the number of calls to FCN
4062 ! with IFLAG = 2.
4063 !
4064 ! Output, integer ( kind = 4 ) IPVT(N), defines a permutation matrix P such
4065 ! that JAC * P = Q * R, where JAC is the final calculated jacobian, Q is
4066 ! orthogonal (not stored), and R is upper triangular.
4067 ! Column J of P is column IPVT(J) of the identity matrix.
4068 !
4069 ! Output, real ( kind = 8 ) QTF(N), contains the first N elements of Q'*FVEC.
4070 !
4071  implicit none
4072 
4073  integer(kind=4) ldfjac
4074  integer(kind=4) m
4075  integer(kind=4) n
4076 
4077  real(kind=8) actred
4078  real(kind=8) delta
4079  real(kind=8) diag(n)
4080  real(kind=8) dirder
4081  real(kind=8) enorm
4082  real(kind=8) epsmch
4083  real(kind=8) factor
4084  external fcn
4085  real(kind=8) fjac(ldfjac, n)
4086  real(kind=8) fnorm
4087  real(kind=8) fnorm1
4088  real(kind=8) ftol
4089  real(kind=8) fvec(m)
4090  real(kind=8) gnorm
4091  real(kind=8) gtol
4092  integer(kind=4) i
4093  integer(kind=4) iflag
4094  integer(kind=4) info
4095  integer(kind=4) ipvt(n)
4096  integer(kind=4) iter
4097  integer(kind=4) j
4098  integer(kind=4) l
4099  integer(kind=4) maxfev
4100  integer(kind=4) mode
4101  integer(kind=4) nfev
4102  integer(kind=4) njev
4103  integer(kind=4) nprint
4104  real(kind=8) par
4105  logical pivot
4106  real(kind=8) pnorm
4107  real(kind=8) prered
4108  real(kind=8) qtf(n)
4109  real(kind=8) ratio
4110  logical sing
4111  real(kind=8) sum2
4112  real(kind=8) temp
4113  real(kind=8) temp1
4114  real(kind=8) temp2
4115  real(kind=8) wa1(n)
4116  real(kind=8) wa2(n)
4117  real(kind=8) wa3(n)
4118  real(kind=8) wa4(m)
4119  real(kind=8) x(n)
4120  real(kind=8) xnorm
4121  real(kind=8) xtol
4122 
4123  epsmch = epsilon(epsmch)
4124 
4125  info = 0
4126  iflag = 0
4127  nfev = 0
4128  njev = 0
4129  xnorm = 0
4130 !
4131 ! Check the input parameters for errors.
4132 !
4133  if (n <= 0) then
4134  go to 340
4135  else if (m < n) then
4136  go to 340
4137  else if (ldfjac < n) then
4138  go to 340
4139  else if (ftol < 0.0d+00) then
4140  go to 340
4141  else if (xtol < 0.0d+00) then
4142  go to 340
4143  else if (gtol < 0.0d+00) then
4144  go to 340
4145  else if (maxfev <= 0) then
4146  go to 340
4147  else if (factor <= 0.0d+00) then
4148  go to 340
4149  end if
4150 
4151  if (mode == 2) then
4152  do j = 1, n
4153  if (diag(j) <= 0.0d+00) then
4154  go to 340
4155  end if
4156  end do
4157  end if
4158 !
4159 ! Evaluate the function at the starting point and calculate its norm.
4160 !
4161  iflag = 1
4162  call fcn(m, n, x, fvec, wa3, iflag)
4163  nfev = 1
4164 
4165  if (iflag < 0) then
4166  go to 340
4167  end if
4168 
4169  fnorm = enorm(m, fvec)
4170 !
4171 ! Initialize Levenberg-Marquardt parameter and iteration counter.
4172 !
4173  par = 0.0d+00
4174  iter = 1
4175 !
4176 ! Beginning of the outer loop.
4177 !
4178 30 continue
4179 !
4180 ! If requested, call FCN to enable printing of iterates.
4181 !
4182  if (0 < nprint) then
4183  iflag = 0
4184  if (mod(iter - 1, nprint) == 0) then
4185  call fcn(m, n, x, fvec, wa3, iflag)
4186  end if
4187  if (iflag < 0) then
4188  go to 340
4189  end if
4190  end if
4191 !
4192 ! Compute the QR factorization of the jacobian matrix calculated one row
4193 ! at a time, while simultaneously forming Q'* FVEC and storing
4194 ! the first N components in QTF.
4195 !
4196  qtf(1:n) = 0.0d+00
4197  fjac(1:n, 1:n) = 0.0d+00
4198  iflag = 2
4199 
4200  do i = 1, m
4201  call fcn(m, n, x, fvec, wa3, iflag)
4202  if (iflag < 0) then
4203  go to 340
4204  end if
4205  temp = fvec(i)
4206  call rwupdt(n, fjac, ldfjac, wa3, qtf, temp, wa1, wa2)
4207  iflag = iflag + 1
4208  end do
4209 
4210  njev = njev + 1
4211 !
4212 ! If the jacobian is rank deficient, call QRFAC to
4213 ! reorder its columns and update the components of QTF.
4214 !
4215  sing = .false.
4216  do j = 1, n
4217  if (fjac(j, j) == 0.0d+00) then
4218  sing = .true.
4219  end if
4220  ipvt(j) = j
4221  wa2(j) = enorm(j, fjac(1, j))
4222  end do
4223 
4224  if (sing) then
4225 
4226  pivot = .true.
4227  call qrfac(n, n, fjac, ldfjac, pivot, ipvt, n, wa1, wa2)
4228 
4229  do j = 1, n
4230 
4231  if (fjac(j, j) /= 0.0d+00) then
4232 
4233  sum2 = dot_product(qtf(j:n), fjac(j:n, j))
4234  temp = -sum2/fjac(j, j)
4235  qtf(j:n) = qtf(j:n) + fjac(j:n, j)*temp
4236 
4237  end if
4238 
4239  fjac(j, j) = wa1(j)
4240 
4241  end do
4242 
4243  end if
4244 !
4245 ! On the first iteration
4246 ! if mode is 1,
4247 ! scale according to the norms of the columns of the initial jacobian.
4248 ! calculate the norm of the scaled X,
4249 ! initialize the step bound delta.
4250 !
4251  if (iter == 1) then
4252 
4253  if (mode /= 2) then
4254 
4255  diag(1:n) = wa2(1:n)
4256  do j = 1, n
4257  if (wa2(j) == 0.0d+00) then
4258  diag(j) = 1.0d+00
4259  end if
4260  end do
4261 
4262  end if
4263 
4264  wa3(1:n) = diag(1:n)*x(1:n)
4265  xnorm = enorm(n, wa3)
4266  delta = factor*xnorm
4267  if (delta == 0.0d+00) then
4268  delta = factor
4269  end if
4270 
4271  end if
4272 !
4273 ! Compute the norm of the scaled gradient.
4274 !
4275  gnorm = 0.0d+00
4276 
4277  if (fnorm /= 0.0d+00) then
4278 
4279  do j = 1, n
4280  l = ipvt(j)
4281  if (wa2(l) /= 0.0d+00) then
4282  sum2 = dot_product(qtf(1:j), fjac(1:j, j))/fnorm
4283  gnorm = max(gnorm, abs(sum2/wa2(l)))
4284  end if
4285  end do
4286 
4287  end if
4288 !
4289 ! Test for convergence of the gradient norm.
4290 !
4291  if (gnorm <= gtol) then
4292  info = 4
4293  go to 340
4294  end if
4295 !
4296 ! Rescale if necessary.
4297 !
4298  if (mode /= 2) then
4299  do j = 1, n
4300  diag(j) = max(diag(j), wa2(j))
4301  end do
4302  end if
4303 !
4304 ! Beginning of the inner loop.
4305 !
4306 240 continue
4307 !
4308 ! Determine the Levenberg-Marquardt parameter.
4309 !
4310  call lmpar(n, fjac, ldfjac, ipvt, diag, qtf, delta, par, wa1, wa2)
4311 !
4312 ! Store the direction P and X + P.
4313 ! Calculate the norm of P.
4314 !
4315  wa1(1:n) = -wa1(1:n)
4316  wa2(1:n) = x(1:n) + wa1(1:n)
4317  wa3(1:n) = diag(1:n)*wa1(1:n)
4318  pnorm = enorm(n, wa3)
4319 !
4320 ! On the first iteration, adjust the initial step bound.
4321 !
4322  if (iter == 1) then
4323  delta = min(delta, pnorm)
4324  end if
4325 !
4326 ! Evaluate the function at X + P and calculate its norm.
4327 !
4328  iflag = 1
4329  call fcn(m, n, wa2, wa4, wa3, iflag)
4330  nfev = nfev + 1
4331  if (iflag < 0) then
4332  go to 340
4333  end if
4334  fnorm1 = enorm(m, wa4)
4335 !
4336 ! Compute the scaled actual reduction.
4337 !
4338  if (0.1d+00*fnorm1 < fnorm) then
4339  actred = 1.0d+00 - (fnorm1/fnorm)**2
4340  else
4341  actred = -1.0d+00
4342  end if
4343 !
4344 ! Compute the scaled predicted reduction and
4345 ! the scaled directional derivative.
4346 !
4347  do j = 1, n
4348  wa3(j) = 0.0d+00
4349  l = ipvt(j)
4350  temp = wa1(l)
4351  wa3(1:j) = wa3(1:j) + fjac(1:j, j)*temp
4352  end do
4353 
4354  temp1 = enorm(n, wa3)/fnorm
4355  temp2 = (sqrt(par)*pnorm)/fnorm
4356  prered = temp1**2 + temp2**2/0.5d+00
4357  dirder = -(temp1**2 + temp2**2)
4358 !
4359 ! Compute the ratio of the actual to the predicted reduction.
4360 !
4361  if (prered /= 0.0d+00) then
4362  ratio = actred/prered
4363  else
4364  ratio = 0.0d+00
4365  end if
4366 !
4367 ! Update the step bound.
4368 !
4369  if (ratio <= 0.25d+00) then
4370 
4371  if (actred >= 0.0d+00) then
4372  temp = 0.5d+00
4373  else
4374  temp = 0.5d+00*dirder/(dirder + 0.5d+00*actred)
4375  end if
4376 
4377  if (0.1d+00*fnorm1 >= fnorm .or. temp < 0.1d+00) then
4378  temp = 0.1d+00
4379  end if
4380 
4381  delta = temp*min(delta, pnorm/0.1d+00)
4382  par = par/temp
4383 
4384  else
4385 
4386  if (par == 0.0d+00 .or. ratio >= 0.75d+00) then
4387  delta = pnorm/0.5d+00
4388  par = 0.5d+00*par
4389  end if
4390 
4391  end if
4392 !
4393 ! Test for successful iteration.
4394 !
4395  if (ratio >= 0.0001d+00) then
4396  x(1:n) = wa2(1:n)
4397  wa2(1:n) = diag(1:n)*x(1:n)
4398  fvec(1:m) = wa4(1:m)
4399  xnorm = enorm(n, wa2)
4400  fnorm = fnorm1
4401  iter = iter + 1
4402  end if
4403 !
4404 ! Tests for convergence, termination and stringent tolerances.
4405 !
4406  if (abs(actred) <= ftol .and. prered <= ftol &
4407  .and. 0.5d+00*ratio <= 1.0d+00) then
4408  info = 1
4409  end if
4410 
4411  if (delta <= xtol*xnorm) then
4412  info = 2
4413  end if
4414 
4415  if (abs(actred) <= ftol .and. prered <= ftol &
4416  .and. 0.5d+00*ratio <= 1.0d+00 .and. info == 2) then
4417  info = 3
4418  end if
4419 
4420  if (info /= 0) then
4421  go to 340
4422  end if
4423 
4424  if (nfev >= maxfev) then
4425  info = 5
4426  end if
4427 
4428  if (abs(actred) <= epsmch .and. prered <= epsmch &
4429  .and. 0.5d+00*ratio <= 1.0d+00) then
4430  info = 6
4431  end if
4432 
4433  if (delta <= epsmch*xnorm) then
4434  info = 7
4435  end if
4436 
4437  if (gnorm <= epsmch) then
4438  info = 8
4439  end if
4440 
4441  if (info /= 0) then
4442  go to 340
4443  end if
4444 !
4445 ! End of the inner loop. Repeat if iteration unsuccessful.
4446 !
4447  if (ratio < 0.0001d+00) then
4448  go to 240
4449  end if
4450 !
4451 ! End of the outer loop.
4452 !
4453  go to 30
4454 
4455 340 continue
4456 !
4457 ! Termination, either normal or user imposed.
4458 !
4459  if (iflag < 0) then
4460  info = iflag
4461  end if
4462 
4463  iflag = 0
4464 
4465  if (0 < nprint) then
4466  call fcn(m, n, x, fvec, wa3, iflag)
4467  end if
4468 
4469  return
subroutine rwupdt(n, r, ldr, w, b, alpha, c, s)
real(kind=8) function enorm(n, x)
subroutine qrfac(m, n, a, lda, pivot, ipvt, lipvt, rdiag, acnorm)
subroutine lmpar(n, r, ldr, ipvt, diag, qtb, delta, par, x, sdiag)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ lmstr1()

subroutine lmstr1 ( external  fcn,
integer(kind=4)  m,
integer(kind=4)  n,
real(kind=8), dimension(n)  x,
real(kind=8), dimension(m)  fvec,
real(kind=8), dimension(ldfjac, n)  fjac,
integer(kind=4)  ldfjac,
real(kind=8)  tol,
integer(kind=4)  info 
)

Definition at line 4472 of file suews_util_minpack.f95.

References lmstr().

4472 
4473 !*****************************************************************************80
4474 !
4475 !! LMSTR1 minimizes M functions in N variables using Levenberg-Marquardt method.
4476 !
4477 ! Discussion:
4478 !
4479 ! LMSTR1 minimizes the sum of the squares of M nonlinear functions in
4480 ! N variables by a modification of the Levenberg-Marquardt algorithm
4481 ! which uses minimal storage.
4482 !
4483 ! This is done by using the more general least-squares solver
4484 ! LMSTR. The user must provide a subroutine which calculates
4485 ! the functions and the rows of the jacobian.
4486 !
4487 ! Licensing:
4488 !
4489 ! This code is distributed under the GNU LGPL license.
4490 !
4491 ! Modified:
4492 !
4493 ! 19 August 2016
4494 !
4495 ! Author:
4496 !
4497 ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom.
4498 ! FORTRAN90 version by John Burkardt.
4499 !
4500 ! Reference:
4501 !
4502 ! Jorge More, Burton Garbow, Kenneth Hillstrom,
4503 ! User Guide for MINPACK-1,
4504 ! Technical Report ANL-80-74,
4505 ! Argonne National Laboratory, 1980.
4506 !
4507 ! Parameters:
4508 !
4509 ! Input, external FCN, the name of the user-supplied subroutine which
4510 ! calculates the functions and the rows of the jacobian.
4511 ! FCN should have the form:
4512 ! subroutine fcn ( m, n, x, fvec, fjrow, iflag )
4513 ! integer ( kind = 4 ) m
4514 ! integer ( kind = 4 ) n
4515 ! real ( kind = 8 ) fjrow(n)
4516 ! real ( kind = 8 ) fvec(m)
4517 ! integer ( kind = 4 ) iflag
4518 ! real ( kind = 8 ) x(n)
4519 ! If IFLAG = 0 on input, then FCN is only being called to allow the user
4520 ! to print out the current iterate.
4521 ! If IFLAG = 1 on input, FCN should calculate the functions at X and
4522 ! return this vector in FVEC.
4523 ! If the input value of IFLAG is I > 1, calculate the (I-1)-st row of
4524 ! the jacobian at X, and return this vector in FJROW.
4525 ! To terminate the algorithm, set the output value of IFLAG negative.
4526 !
4527 ! Input, integer ( kind = 4 ) M, the number of functions.
4528 !
4529 ! Input, integer ( kind = 4 ) N, the number of variables.
4530 ! N must not exceed M.
4531 !
4532 ! Input/output, real ( kind = 8 ) X(N). On input, X must contain an initial
4533 ! estimate of the solution vector. On output X contains the final
4534 ! estimate of the solution vector.
4535 !
4536 ! Output, real ( kind = 8 ) FVEC(M), the functions evaluated at the output X.
4537 !
4538 ! Output, real ( kind = 8 ) FJAC(LDFJAC,N), an N by N array. The upper
4539 ! triangle contains an upper triangular matrix R such that
4540 !
4541 ! P' * ( JAC' * JAC ) * P = R' * R,
4542 !
4543 ! where P is a permutation matrix and JAC is the final calculated
4544 ! jacobian. Column J of P is column IPVT(J) of the identity matrix.
4545 ! The lower triangular part of FJAC contains information generated
4546 ! during the computation of R.
4547 !
4548 ! Input, integer ( kind = 4 ) LDFJAC, the leading dimension of FJAC.
4549 ! LDFJAC must be at least N.
4550 !
4551 ! Input, real ( kind = 8 ) TOL. Termination occurs when the algorithm
4552 ! estimates either that the relative error in the sum of squares is at
4553 ! most TOL or that the relative error between X and the solution is at
4554 ! most TOL. TOL should be nonnegative.
4555 !
4556 ! Output, integer ( kind = 4 ) INFO, error flag. If the user has terminated
4557 ! execution, INFO is set to the (negative) value of IFLAG. See description
4558 ! of FCN. Otherwise, INFO is set as follows:
4559 ! 0, improper input parameters.
4560 ! 1, algorithm estimates that the relative error in the sum of squares
4561 ! is at most TOL.
4562 ! 2, algorithm estimates that the relative error between X and the
4563 ! solution is at most TOL.
4564 ! 3, conditions for INFO = 1 and INFO = 2 both hold.
4565 ! 4, FVEC is orthogonal to the columns of the jacobian to machine precision.
4566 ! 5, number of calls to FCN with IFLAG = 1 has reached 100*(N+1).
4567 ! 6, TOL is too small. No further reduction in the sum of squares
4568 ! is possible.
4569 ! 7, TOL is too small. No further improvement in the approximate
4570 ! solution X is possible.
4571 !
4572  implicit none
4573 
4574  integer(kind=4) ldfjac
4575  integer(kind=4) m
4576  integer(kind=4) n
4577 
4578  real(kind=8) diag(n)
4579  real(kind=8) factor
4580  external fcn
4581  real(kind=8) fjac(ldfjac, n)
4582  real(kind=8) ftol
4583  real(kind=8) fvec(m)
4584  real(kind=8) gtol
4585  integer(kind=4) info
4586  integer(kind=4) ipvt(n)
4587  integer(kind=4) maxfev
4588  integer(kind=4) mode
4589  integer(kind=4) nfev
4590  integer(kind=4) njev
4591  integer(kind=4) nprint
4592  real(kind=8) qtf(n)
4593  real(kind=8) tol
4594  real(kind=8) x(n)
4595  real(kind=8) xtol
4596 
4597  if (n <= 0) then
4598  info = 0
4599  return
4600  end if
4601 
4602  if (m < n) then
4603  info = 0
4604  return
4605  end if
4606 
4607  if (ldfjac < n) then
4608  info = 0
4609  return
4610  end if
4611 
4612  if (tol < 0.0d+00) then
4613  info = 0
4614  return
4615  end if
4616 
4617  fvec(1:n) = 0.0d+00
4618  fjac(1:ldfjac, 1:n) = 0.0d+00
4619  ftol = tol
4620  xtol = tol
4621  gtol = 0.0d+00
4622  maxfev = 100*(n + 1)
4623  diag(1:n) = 0.0d+00
4624  mode = 1
4625  factor = 100.0d+00
4626  nprint = 0
4627  info = 0
4628  nfev = 0
4629  njev = 0
4630  ipvt(1:n) = 0
4631  qtf(1:n) = 0.0d+00
4632 
4633  call lmstr(fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, &
4634  diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf)
4635 
4636  if (info == 8) then
4637  info = 4
4638  end if
4639 
4640  return
subroutine lmstr(fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf)
Here is the call graph for this function:

◆ qform()

subroutine qform ( integer(kind=4)  m,
integer(kind=4)  n,
real(kind=8), dimension(ldq, m)  q,
integer(kind=4)  ldq 
)

Definition at line 4643 of file suews_util_minpack.f95.

Referenced by hybrd(), and hybrj().

4643 
4644 !*****************************************************************************80
4645 !
4646 !! QFORM produces the explicit QR factorization of a matrix.
4647 !
4648 ! Discussion:
4649 !
4650 ! The QR factorization of a matrix is usually accumulated in implicit
4651 ! form, that is, as a series of orthogonal transformations of the
4652 ! original matrix. This routine carries out those transformations,
4653 ! to explicitly exhibit the factorization constructed by QRFAC.
4654 !
4655 ! Licensing:
4656 !
4657 ! This code is distributed under the GNU LGPL license.
4658 !
4659 ! Modified:
4660 !
4661 ! 06 April 2010
4662 !
4663 ! Author:
4664 !
4665 ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom.
4666 ! FORTRAN90 version by John Burkardt.
4667 !
4668 ! Reference:
4669 !
4670 ! Jorge More, Burton Garbow, Kenneth Hillstrom,
4671 ! User Guide for MINPACK-1,
4672 ! Technical Report ANL-80-74,
4673 ! Argonne National Laboratory, 1980.
4674 !
4675 ! Parameters:
4676 !
4677 ! Input, integer ( kind = 4 ) M, is a positive integer input variable set
4678 ! to the number of rows of A and the order of Q.
4679 !
4680 ! Input, integer ( kind = 4 ) N, is a positive integer input variable set
4681 ! to the number of columns of A.
4682 !
4683 ! Input/output, real ( kind = 8 ) Q(LDQ,M). Q is an M by M array.
4684 ! On input the full lower trapezoid in the first min(M,N) columns of Q
4685 ! contains the factored form.
4686 ! On output, Q has been accumulated into a square matrix.
4687 !
4688 ! Input, integer ( kind = 4 ) LDQ, is a positive integer input variable
4689 ! not less than M which specifies the leading dimension of the array Q.
4690 !
4691  implicit none
4692 
4693  integer(kind=4) ldq
4694  integer(kind=4) m
4695  integer(kind=4) n
4696 
4697  integer(kind=4) j
4698  integer(kind=4) k
4699  integer(kind=4) l
4700  integer(kind=4) minmn
4701  real(kind=8) q(ldq, m)
4702  real(kind=8) temp
4703  real(kind=8) wa(m)
4704 
4705  minmn = min(m, n)
4706 
4707  do j = 2, minmn
4708  q(1:j - 1, j) = 0.0d+00
4709  end do
4710 !
4711 ! Initialize remaining columns to those of the identity matrix.
4712 !
4713  q(1:m, n + 1:m) = 0.0d+00
4714 
4715  do j = n + 1, m
4716  q(j, j) = 1.0d+00
4717  end do
4718 !
4719 ! Accumulate Q from its factored form.
4720 !
4721  do l = 1, minmn
4722 
4723  k = minmn - l + 1
4724 
4725  wa(k:m) = q(k:m, k)
4726 
4727  q(k:m, k) = 0.0d+00
4728  q(k, k) = 1.0d+00
4729 
4730  if (wa(k) /= 0.0d+00) then
4731 
4732  do j = k, m
4733  temp = dot_product(wa(k:m), q(k:m, j))/wa(k)
4734  q(k:m, j) = q(k:m, j) - temp*wa(k:m)
4735  end do
4736 
4737  end if
4738 
4739  end do
4740 
4741  return
Here is the caller graph for this function:

◆ qrfac()

subroutine qrfac ( integer(kind=4)  m,
integer(kind=4)  n,
real(kind=8), dimension(lda, n)  a,
integer(kind=4)  lda,
logical  pivot,
integer(kind=4), dimension(lipvt)  ipvt,
integer(kind=4)  lipvt,
real(kind=8), dimension(n)  rdiag,
real(kind=8), dimension(n)  acnorm 
)

Definition at line 4744 of file suews_util_minpack.f95.

Referenced by hybrd(), hybrj(), lmder(), lmdif(), and lmstr().

4744 
4745 !*****************************************************************************80
4746 !
4747 !! QRFAC computes a QR factorization using Householder transformations.
4748 !
4749 ! Discussion:
4750 !
4751 ! This subroutine uses Householder transformations with column
4752 ! pivoting (optional) to compute a QR factorization of the
4753 ! M by N matrix A. That is, QRFAC determines an orthogonal
4754 ! matrix Q, a permutation matrix P, and an upper trapezoidal
4755 ! matrix R with diagonal elements of nonincreasing magnitude,
4756 ! such that A*P = Q*R. The Householder transformation for
4757 ! column K, K = 1,2,...,min(M,N), is of the form
4758 !
4759 ! I - ( 1 / U(K) ) * U * U'
4760 !
4761 ! where U has zeros in the first K-1 positions. The form of
4762 ! this transformation and the method of pivoting first
4763 ! appeared in the corresponding LINPACK subroutine.
4764 !
4765 ! Licensing:
4766 !
4767 ! This code is distributed under the GNU LGPL license.
4768 !
4769 ! Modified:
4770 !
4771 ! 06 April 2010
4772 !
4773 ! Author:
4774 !
4775 ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom.
4776 ! FORTRAN90 version by John Burkardt.
4777 !
4778 ! Reference:
4779 !
4780 ! Jorge More, Burton Garbow, Kenneth Hillstrom,
4781 ! User Guide for MINPACK-1,
4782 ! Technical Report ANL-80-74,
4783 ! Argonne National Laboratory, 1980.
4784 !
4785 ! Parameters:
4786 !
4787 ! Input, integer ( kind = 4 ) M, the number of rows of A.
4788 !
4789 ! Input, integer ( kind = 4 ) N, the number of columns of A.
4790 !
4791 ! Input/output, real ( kind = 8 ) A(LDA,N), the M by N array.
4792 ! On input, A contains the matrix for which the QR factorization is to
4793 ! be computed. On output, the strict upper trapezoidal part of A contains
4794 ! the strict upper trapezoidal part of R, and the lower trapezoidal
4795 ! part of A contains a factored form of Q (the non-trivial elements of
4796 ! the U vectors described above).
4797 !
4798 ! Input, integer ( kind = 4 ) LDA, the leading dimension of A, which must
4799 ! be no less than M.
4800 !
4801 ! Input, logical PIVOT, is TRUE if column pivoting is to be carried out.
4802 !
4803 ! Output, integer ( kind = 4 ) IPVT(LIPVT), defines the permutation matrix P
4804 ! such that A*P = Q*R. Column J of P is column IPVT(J) of the identity
4805 ! matrix. If PIVOT is false, IPVT is not referenced.
4806 !
4807 ! Input, integer ( kind = 4 ) LIPVT, the dimension of IPVT, which should
4808 ! be N if pivoting is used.
4809 !
4810 ! Output, real ( kind = 8 ) RDIAG(N), contains the diagonal elements of R.
4811 !
4812 ! Output, real ( kind = 8 ) ACNORM(N), the norms of the corresponding
4813 ! columns of the input matrix A. If this information is not needed,
4814 ! then ACNORM can coincide with RDIAG.
4815 !
4816  implicit none
4817 
4818  integer(kind=4) lda
4819  integer(kind=4) lipvt
4820  integer(kind=4) m
4821  integer(kind=4) n
4822 
4823  real(kind=8) a(lda, n)
4824  real(kind=8) acnorm(n)
4825  real(kind=8) ajnorm
4826  real(kind=8) enorm
4827  real(kind=8) epsmch
4828  ! integer ( kind = 4 ) i
4829  integer(kind=4) i4_temp
4830  integer(kind=4) ipvt(lipvt)
4831  integer(kind=4) j
4832  integer(kind=4) k
4833  integer(kind=4) kmax
4834  integer(kind=4) minmn
4835  logical pivot
4836  real(kind=8) r8_temp(m)
4837  real(kind=8) rdiag(n)
4838  real(kind=8) temp
4839  real(kind=8) wa(n)
4840 
4841  epsmch = epsilon(epsmch)
4842 !
4843 ! Compute the initial column norms and initialize several arrays.
4844 !
4845  do j = 1, n
4846  acnorm(j) = enorm(m, a(1:m, j))
4847  end do
4848 
4849  rdiag(1:n) = acnorm(1:n)
4850  wa(1:n) = acnorm(1:n)
4851 
4852  if (pivot) then
4853  do j = 1, n
4854  ipvt(j) = j
4855  end do
4856  end if
4857 !
4858 ! Reduce A to R with Householder transformations.
4859 !
4860  minmn = min(m, n)
4861 
4862  do j = 1, minmn
4863 !
4864 ! Bring the column of largest norm into the pivot position.
4865 !
4866  if (pivot) then
4867 
4868  kmax = j
4869 
4870  do k = j, n
4871  if (rdiag(kmax) < rdiag(k)) then
4872  kmax = k
4873  end if
4874  end do
4875 
4876  if (kmax /= j) then
4877 
4878  r8_temp(1:m) = a(1:m, j)
4879  a(1:m, j) = a(1:m, kmax)
4880  a(1:m, kmax) = r8_temp(1:m)
4881 
4882  rdiag(kmax) = rdiag(j)
4883  wa(kmax) = wa(j)
4884 
4885  i4_temp = ipvt(j)
4886  ipvt(j) = ipvt(kmax)
4887  ipvt(kmax) = i4_temp
4888 
4889  end if
4890 
4891  end if
4892 !
4893 ! Compute the Householder transformation to reduce the
4894 ! J-th column of A to a multiple of the J-th unit vector.
4895 !
4896  ajnorm = enorm(m - j + 1, a(j, j))
4897 
4898  if (ajnorm /= 0.0d+00) then
4899 
4900  if (a(j, j) < 0.0d+00) then
4901  ajnorm = -ajnorm
4902  end if
4903 
4904  a(j:m, j) = a(j:m, j)/ajnorm
4905  a(j, j) = a(j, j) + 1.0d+00
4906 !
4907 ! Apply the transformation to the remaining columns and update the norms.
4908 !
4909  do k = j + 1, n
4910 
4911  temp = dot_product(a(j:m, j), a(j:m, k))/a(j, j)
4912 
4913  a(j:m, k) = a(j:m, k) - temp*a(j:m, j)
4914 
4915  if (pivot .and. rdiag(k) /= 0.0d+00) then
4916 
4917  temp = a(j, k)/rdiag(k)
4918  rdiag(k) = rdiag(k)*sqrt(max(0.0d+00, 1.0d+00 - temp**2))
4919 
4920  if (0.05d+00*(rdiag(k)/wa(k))**2 <= epsmch) then
4921  rdiag(k) = enorm(m - j, a(j + 1, k))
4922  wa(k) = rdiag(k)
4923  end if
4924 
4925  end if
4926 
4927  end do
4928 
4929  end if
4930 
4931  rdiag(j) = -ajnorm
4932 
4933  end do
4934 
4935  return
real(kind(1d0)) kmax
real(kind=8) function enorm(n, x)
Here is the caller graph for this function:

◆ qrsolv()

subroutine qrsolv ( integer(kind=4)  n,
real(kind=8), dimension(ldr, n)  r,
integer(kind=4)  ldr,
integer(kind=4), dimension(n)  ipvt,
real(kind=8), dimension(n)  diag,
real(kind=8), dimension(n)  qtb,
real(kind=8), dimension(n)  x,
real(kind=8), dimension(n)  sdiag 
)

Definition at line 4938 of file suews_util_minpack.f95.

Referenced by lmpar().

4938 
4939 !*****************************************************************************80
4940 !
4941 !! QRSOLV solves a rectangular linear system A*x=b in the least squares sense.
4942 !
4943 ! Discussion:
4944 !
4945 ! Given an M by N matrix A, an N by N diagonal matrix D,
4946 ! and an M-vector B, the problem is to determine an X which
4947 ! solves the system
4948 !
4949 ! A*X = B
4950 ! D*X = 0
4951 !
4952 ! in the least squares sense.
4953 !
4954 ! This subroutine completes the solution of the problem
4955 ! if it is provided with the necessary information from the
4956 ! QR factorization, with column pivoting, of A. That is, if
4957 ! Q*P = Q*R, where P is a permutation matrix, Q has orthogonal
4958 ! columns, and R is an upper triangular matrix with diagonal
4959 ! elements of nonincreasing magnitude, then QRSOLV expects
4960 ! the full upper triangle of R, the permutation matrix p,
4961 ! and the first N components of Q'*B.
4962 !
4963 ! The system is then equivalent to
4964 !
4965 ! R*Z = Q'*B
4966 ! P'*D*P*Z = 0
4967 !
4968 ! where X = P*Z. If this system does not have full rank,
4969 ! then a least squares solution is obtained. On output QRSOLV
4970 ! also provides an upper triangular matrix S such that
4971 !
4972 ! P'*(A'*A + D*D)*P = S'*S.
4973 !
4974 ! S is computed within QRSOLV and may be of separate interest.
4975 !
4976 ! Licensing:
4977 !
4978 ! This code is distributed under the GNU LGPL license.
4979 !
4980 ! Modified:
4981 !
4982 ! 06 April 2010
4983 !
4984 ! Author:
4985 !
4986 ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom.
4987 ! FORTRAN90 version by John Burkardt.
4988 !
4989 ! Reference:
4990 !
4991 ! Jorge More, Burton Garbow, Kenneth Hillstrom,
4992 ! User Guide for MINPACK-1,
4993 ! Technical Report ANL-80-74,
4994 ! Argonne National Laboratory, 1980.
4995 !
4996 ! Parameters:
4997 !
4998 ! Input, integer ( kind = 4 ) N, the order of R.
4999 !
5000 ! Input/output, real ( kind = 8 ) R(LDR,N), the N by N matrix.
5001 ! On input the full upper triangle must contain the full upper triangle
5002 ! of the matrix R. On output the full upper triangle is unaltered, and
5003 ! the strict lower triangle contains the strict upper triangle
5004 ! (transposed) of the upper triangular matrix S.
5005 !
5006 ! Input, integer ( kind = 4 ) LDR, the leading dimension of R, which must be
5007 ! at least N.
5008 !
5009 ! Input, integer ( kind = 4 ) IPVT(N), defines the permutation matrix P such
5010 ! that A*P = Q*R. Column J of P is column IPVT(J) of the identity matrix.
5011 !
5012 ! Input, real ( kind = 8 ) DIAG(N), the diagonal elements of the matrix D.
5013 !
5014 ! Input, real ( kind = 8 ) QTB(N), the first N elements of the vector Q'*B.
5015 !
5016 ! Output, real ( kind = 8 ) X(N), the least squares solution.
5017 !
5018 ! Output, real ( kind = 8 ) SDIAG(N), the diagonal elements of the upper
5019 ! triangular matrix S.
5020 !
5021  implicit none
5022 
5023  integer(kind=4) ldr
5024  integer(kind=4) n
5025 
5026  real(kind=8) c
5027  real(kind=8) cotan
5028  real(kind=8) diag(n)
5029  integer(kind=4) i
5030  integer(kind=4) ipvt(n)
5031  integer(kind=4) j
5032  integer(kind=4) k
5033  integer(kind=4) l
5034  integer(kind=4) nsing
5035  real(kind=8) qtb(n)
5036  real(kind=8) qtbpj
5037  real(kind=8) r(ldr, n)
5038  real(kind=8) s
5039  real(kind=8) sdiag(n)
5040  real(kind=8) sum2
5041  real(kind=8) t
5042  real(kind=8) temp
5043  real(kind=8) wa(n)
5044  real(kind=8) x(n)
5045 !
5046 ! Copy R and Q'*B to preserve input and initialize S.
5047 !
5048 ! In particular, save the diagonal elements of R in X.
5049 !
5050  do j = 1, n
5051  r(j:n, j) = r(j, j:n)
5052  x(j) = r(j, j)
5053  end do
5054 
5055  wa(1:n) = qtb(1:n)
5056 !
5057 ! Eliminate the diagonal matrix D using a Givens rotation.
5058 !
5059  do j = 1, n
5060 !
5061 ! Prepare the row of D to be eliminated, locating the
5062 ! diagonal element using P from the QR factorization.
5063 !
5064  l = ipvt(j)
5065 
5066  if (diag(l) /= 0.0d+00) then
5067 
5068  sdiag(j:n) = 0.0d+00
5069  sdiag(j) = diag(l)
5070 !
5071 ! The transformations to eliminate the row of D
5072 ! modify only a single element of Q'*B
5073 ! beyond the first N, which is initially zero.
5074 !
5075  qtbpj = 0.0d+00
5076 
5077  do k = j, n
5078 !
5079 ! Determine a Givens rotation which eliminates the
5080 ! appropriate element in the current row of D.
5081 !
5082  if (sdiag(k) /= 0.0d+00) then
5083 
5084  if (abs(r(k, k)) < abs(sdiag(k))) then
5085  cotan = r(k, k)/sdiag(k)
5086  s = 0.5d+00/sqrt(0.25d+00 + 0.25d+00*cotan**2)
5087  c = s*cotan
5088  else
5089  t = sdiag(k)/r(k, k)
5090  c = 0.5d+00/sqrt(0.25d+00 + 0.25d+00*t**2)
5091  s = c*t
5092  end if
5093 !
5094 ! Compute the modified diagonal element of R and
5095 ! the modified element of (Q'*B,0).
5096 !
5097  r(k, k) = c*r(k, k) + s*sdiag(k)
5098  temp = c*wa(k) + s*qtbpj
5099  qtbpj = -s*wa(k) + c*qtbpj
5100  wa(k) = temp
5101 !
5102 ! Accumulate the tranformation in the row of S.
5103 !
5104  do i = k + 1, n
5105  temp = c*r(i, k) + s*sdiag(i)
5106  sdiag(i) = -s*r(i, k) + c*sdiag(i)
5107  r(i, k) = temp
5108  end do
5109 
5110  end if
5111 
5112  end do
5113 
5114  end if
5115 !
5116 ! Store the diagonal element of S and restore
5117 ! the corresponding diagonal element of R.
5118 !
5119  sdiag(j) = r(j, j)
5120  r(j, j) = x(j)
5121 
5122  end do
5123 !
5124 ! Solve the triangular system for Z. If the system is
5125 ! singular, then obtain a least squares solution.
5126 !
5127  nsing = n
5128 
5129  do j = 1, n
5130 
5131  if (sdiag(j) == 0.0d+00 .and. nsing == n) then
5132  nsing = j - 1
5133  end if
5134 
5135  if (nsing < n) then
5136  wa(j) = 0.0d+00
5137  end if
5138 
5139  end do
5140 
5141  do j = nsing, 1, -1
5142  sum2 = dot_product(wa(j + 1:nsing), r(j + 1:nsing, j))
5143  wa(j) = (wa(j) - sum2)/sdiag(j)
5144  end do
5145 !
5146 ! Permute the components of Z back to components of X.
5147 !
5148  do j = 1, n
5149  l = ipvt(j)
5150  x(l) = wa(j)
5151  end do
5152 
5153  return
Here is the caller graph for this function:

◆ r1mpyq()

subroutine r1mpyq ( integer(kind=4)  m,
integer(kind=4)  n,
real(kind=8), dimension(lda, n)  a,
integer(kind=4)  lda,
real(kind=8), dimension(n)  v,
real(kind=8), dimension(n)  w 
)

Definition at line 5156 of file suews_util_minpack.f95.

Referenced by hybrd(), and hybrj().

5156 
5157 !*****************************************************************************80
5158 !
5159 !! R1MPYQ computes A*Q, where Q is the product of Householder transformations.
5160 !
5161 ! Discussion:
5162 !
5163 ! Given an M by N matrix A, this subroutine computes A*Q where
5164 ! Q is the product of 2*(N - 1) transformations
5165 !
5166 ! GV(N-1)*...*GV(1)*GW(1)*...*GW(N-1)
5167 !
5168 ! and GV(I), GW(I) are Givens rotations in the (I,N) plane which
5169 ! eliminate elements in the I-th and N-th planes, respectively.
5170 ! Q itself is not given, rather the information to recover the
5171 ! GV, GW rotations is supplied.
5172 !
5173 ! Licensing:
5174 !
5175 ! This code is distributed under the GNU LGPL license.
5176 !
5177 ! Modified:
5178 !
5179 ! 06 April 2010
5180 !
5181 ! Author:
5182 !
5183 ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom.
5184 ! FORTRAN90 version by John Burkardt.
5185 !
5186 ! Reference:
5187 !
5188 ! Jorge More, Burton Garbow, Kenneth Hillstrom,
5189 ! User Guide for MINPACK-1,
5190 ! Technical Report ANL-80-74,
5191 ! Argonne National Laboratory, 1980.
5192 !
5193 ! Parameters:
5194 !
5195 ! Input, integer ( kind = 4 ) M, the number of rows of A.
5196 !
5197 ! Input, integer ( kind = 4 ) N, the number of columns of A.
5198 !
5199 ! Input/output, real ( kind = 8 ) A(LDA,N), the M by N array.
5200 ! On input, the matrix A to be postmultiplied by the orthogonal matrix Q.
5201 ! On output, the value of A*Q.
5202 !
5203 ! Input, integer ( kind = 4 ) LDA, the leading dimension of A, which must not
5204 ! be less than M.
5205 !
5206 ! Input, real ( kind = 8 ) V(N), W(N), contain the information necessary
5207 ! to recover the Givens rotations GV and GW.
5208 !
5209  implicit none
5210 
5211  integer(kind=4) lda
5212  integer(kind=4) m
5213  integer(kind=4) n
5214 
5215  real(kind=8) a(lda, n)
5216  real(kind=8) c
5217  integer(kind=4) i
5218  integer(kind=4) j
5219  real(kind=8) s
5220  real(kind=8) temp
5221  real(kind=8) v(n)
5222  real(kind=8) w(n)
5223 !
5224 ! Apply the first set of Givens rotations to A.
5225 !
5226  do j = n - 1, 1, -1
5227 
5228  if (1.0d+00 < abs(v(j))) then
5229  c = 1.0d+00/v(j)
5230  s = sqrt(1.0d+00 - c**2)
5231  else
5232  s = v(j)
5233  c = sqrt(1.0d+00 - s**2)
5234  end if
5235 
5236  do i = 1, m
5237  temp = c*a(i, j) - s*a(i, n)
5238  a(i, n) = s*a(i, j) + c*a(i, n)
5239  a(i, j) = temp
5240  end do
5241 
5242  end do
5243 !
5244 ! Apply the second set of Givens rotations to A.
5245 !
5246  do j = 1, n - 1
5247 
5248  if (abs(w(j)) > 1.0d+00) then
5249  c = 1.0d+00/w(j)
5250  s = sqrt(1.0d+00 - c**2)
5251  else
5252  s = w(j)
5253  c = sqrt(1.0d+00 - s**2)
5254  end if
5255 
5256  do i = 1, m
5257  temp = c*a(i, j) + s*a(i, n)
5258  a(i, n) = -s*a(i, j) + c*a(i, n)
5259  a(i, j) = temp
5260  end do
5261 
5262  end do
5263 
5264  return
Here is the caller graph for this function:

◆ r1updt()

subroutine r1updt ( integer(kind=4)  m,
integer(kind=4)  n,
real(kind=8), dimension(ls)  s,
integer(kind=4)  ls,
real(kind=8), dimension(m)  u,
real(kind=8), dimension(n)  v,
real(kind=8), dimension(m)  w,
logical  sing 
)

Definition at line 5267 of file suews_util_minpack.f95.

Referenced by hybrd(), and hybrj().

5267 
5268 !*****************************************************************************80
5269 !
5270 !! R1UPDT re-triangularizes a matrix after a rank one update.
5271 !
5272 ! Discussion:
5273 !
5274 ! Given an M by N lower trapezoidal matrix S, an M-vector U, and an
5275 ! N-vector V, the problem is to determine an orthogonal matrix Q such that
5276 !
5277 ! (S + U * V' ) * Q
5278 !
5279 ! is again lower trapezoidal.
5280 !
5281 ! This subroutine determines Q as the product of 2 * (N - 1)
5282 ! transformations
5283 !
5284 ! GV(N-1)*...*GV(1)*GW(1)*...*GW(N-1)
5285 !
5286 ! where GV(I), GW(I) are Givens rotations in the (I,N) plane
5287 ! which eliminate elements in the I-th and N-th planes,
5288 ! respectively. Q itself is not accumulated, rather the
5289 ! information to recover the GV and GW rotations is returned.
5290 !
5291 ! Licensing:
5292 !
5293 ! This code is distributed under the GNU LGPL license.
5294 !
5295 ! Modified:
5296 !
5297 ! 06 April 2010
5298 !
5299 ! Author:
5300 !
5301 ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom.
5302 ! FORTRAN90 version by John Burkardt.
5303 !
5304 ! Reference:
5305 !
5306 ! Jorge More, Burton Garbow, Kenneth Hillstrom,
5307 ! User Guide for MINPACK-1,
5308 ! Technical Report ANL-80-74,
5309 ! Argonne National Laboratory, 1980.
5310 !
5311 ! Parameters:
5312 !
5313 ! Input, integer ( kind = 4 ) M, the number of rows of S.
5314 !
5315 ! Input, integer ( kind = 4 ) N, the number of columns of S.
5316 ! N must not exceed M.
5317 !
5318 ! Input/output, real ( kind = 8 ) S(LS). On input, the lower trapezoidal
5319 ! matrix S stored by columns. On output S contains the lower trapezoidal
5320 ! matrix produced as described above.
5321 !
5322 ! Input, integer ( kind = 4 ) LS, the length of the S array. LS must be at
5323 ! least (N*(2*M-N+1))/2.
5324 !
5325 ! Input, real ( kind = 8 ) U(M), the U vector.
5326 !
5327 ! Input/output, real ( kind = 8 ) V(N). On input, V must contain the
5328 ! vector V. On output V contains the information necessary to recover the
5329 ! Givens rotations GV described above.
5330 !
5331 ! Output, real ( kind = 8 ) W(M), contains information necessary to
5332 ! recover the Givens rotations GW described above.
5333 !
5334 ! Output, logical SING, is set to TRUE if any of the diagonal elements
5335 ! of the output S are zero. Otherwise SING is set FALSE.
5336 !
5337  implicit none
5338 
5339  integer(kind=4) ls
5340  integer(kind=4) m
5341  integer(kind=4) n
5342 
5343  real(kind=8) cos
5344  real(kind=8) cotan
5345  real(kind=8) giant
5346  integer(kind=4) i
5347  integer(kind=4) j
5348  integer(kind=4) jj
5349  integer(kind=4) l
5350  real(kind=8) s(ls)
5351  real(kind=8) sin
5352  logical sing
5353  real(kind=8) tan
5354  real(kind=8) tau
5355  real(kind=8) temp
5356  real(kind=8) u(m)
5357  real(kind=8) v(n)
5358  real(kind=8) w(m)
5359 !
5360 ! GIANT is the largest magnitude.
5361 !
5362  giant = huge(giant)
5363 !
5364 ! Initialize the diagonal element pointer.
5365 !
5366  jj = (n*(2*m - n + 1))/2 - (m - n)
5367 !
5368 ! Move the nontrivial part of the last column of S into W.
5369 !
5370  l = jj
5371  do i = n, m
5372  w(i) = s(l)
5373  l = l + 1
5374  end do
5375 !
5376 ! Rotate the vector V into a multiple of the N-th unit vector
5377 ! in such a way that a spike is introduced into W.
5378 !
5379  do j = n - 1, 1, -1
5380 
5381  jj = jj - (m - j + 1)
5382  w(j) = 0.0d+00
5383 
5384  if (v(j) /= 0.0d+00) then
5385 !
5386 ! Determine a Givens rotation which eliminates the
5387 ! J-th element of V.
5388 !
5389  if (abs(v(n)) < abs(v(j))) then
5390  cotan = v(n)/v(j)
5391  sin = 0.5d+00/sqrt(0.25d+00 + 0.25d+00*cotan**2)
5392  cos = sin*cotan
5393  tau = 1.0d+00
5394  if (abs(cos)*giant > 1.0d+00) then
5395  tau = 1.0d+00/cos
5396  end if
5397  else
5398  tan = v(j)/v(n)
5399  cos = 0.5d+00/sqrt(0.25d+00 + 0.25d+00*tan**2)
5400  sin = cos*tan
5401  tau = sin
5402  end if
5403 !
5404 ! Apply the transformation to V and store the information
5405 ! necessary to recover the Givens rotation.
5406 !
5407  v(n) = sin*v(j) + cos*v(n)
5408  v(j) = tau
5409 !
5410 ! Apply the transformation to S and extend the spike in W.
5411 !
5412  l = jj
5413  do i = j, m
5414  temp = cos*s(l) - sin*w(i)
5415  w(i) = sin*s(l) + cos*w(i)
5416  s(l) = temp
5417  l = l + 1
5418  end do
5419 
5420  end if
5421 
5422  end do
5423 !
5424 ! Add the spike from the rank 1 update to W.
5425 !
5426  w(1:m) = w(1:m) + v(n)*u(1:m)
5427 !
5428 ! Eliminate the spike.
5429 !
5430  sing = .false.
5431 
5432  do j = 1, n - 1
5433 
5434  if (w(j) /= 0.0d+00) then
5435 !
5436 ! Determine a Givens rotation which eliminates the
5437 ! J-th element of the spike.
5438 !
5439  if (abs(s(jj)) < abs(w(j))) then
5440 
5441  cotan = s(jj)/w(j)
5442  sin = 0.5d+00/sqrt(0.25d+00 + 0.25d+00*cotan**2)
5443  cos = sin*cotan
5444 
5445  if (1.0d+00 < abs(cos)*giant) then
5446  tau = 1.0d+00/cos
5447  else
5448  tau = 1.0d+00
5449  end if
5450 
5451  else
5452 
5453  tan = w(j)/s(jj)
5454  cos = 0.5d+00/sqrt(0.25d+00 + 0.25d+00*tan**2)
5455  sin = cos*tan
5456  tau = sin
5457 
5458  end if
5459 !
5460 ! Apply the transformation to S and reduce the spike in W.
5461 !
5462  l = jj
5463  do i = j, m
5464  temp = cos*s(l) + sin*w(i)
5465  w(i) = -sin*s(l) + cos*w(i)
5466  s(l) = temp
5467  l = l + 1
5468  end do
5469 !
5470 ! Store the information necessary to recover the Givens rotation.
5471 !
5472  w(j) = tau
5473 
5474  end if
5475 !
5476 ! Test for zero diagonal elements in the output S.
5477 !
5478  if (s(jj) == 0.0d+00) then
5479  sing = .true.
5480  end if
5481 
5482  jj = jj + (m - j + 1)
5483 
5484  end do
5485 !
5486 ! Move W back into the last column of the output S.
5487 !
5488  l = jj
5489  do i = n, m
5490  s(l) = w(i)
5491  l = l + 1
5492  end do
5493 
5494  if (s(jj) == 0.0d+00) then
5495  sing = .true.
5496  end if
5497 
5498  return
Here is the caller graph for this function:

◆ r8vec_print()

subroutine r8vec_print ( integer(kind=4)  n,
real(kind=8), dimension(n)  a,
character(len=*)  title 
)

Definition at line 5501 of file suews_util_minpack.f95.

Referenced by anohm_module::anohm_fccal(), and anohm_module::fcnbo().

5501 
5502 !*****************************************************************************80
5503 !
5504 !! R8VEC_PRINT prints an R8VEC.
5505 !
5506 ! Discussion:
5507 !
5508 ! An R8VEC is a vector of R8's.
5509 !
5510 ! Licensing:
5511 !
5512 ! This code is distributed under the GNU LGPL license.
5513 !
5514 ! Modified:
5515 !
5516 ! 22 August 2000
5517 !
5518 ! Author:
5519 !
5520 ! John Burkardt
5521 !
5522 ! Parameters:
5523 !
5524 ! Input, integer ( kind = 4 ) N, the number of components of the vector.
5525 !
5526 ! Input, real ( kind = 8 ) A(N), the vector to be printed.
5527 !
5528 ! Input, character ( len = * ) TITLE, a title.
5529 !
5530  implicit none
5531 
5532  integer(kind=4) n
5533 
5534  real(kind=8) a(n)
5535  integer(kind=4) i
5536  character(len=*) title
5537 
5538  write (*, '(a)') ' '
5539  write (*, '(a)') trim(title)
5540  write (*, '(a)') ' '
5541  do i = 1, n
5542  write (*, '(2x,i8,2x,g16.8)') i, a(i)
5543  end do
5544 
5545  return
Here is the caller graph for this function:

◆ rwupdt()

subroutine rwupdt ( integer(kind=4)  n,
real(kind=8), dimension(ldr, n)  r,
integer(kind=4)  ldr,
real(kind=8), dimension(n)  w,
real(kind=8), dimension(n)  b,
real(kind=8)  alpha,
real(kind=8), dimension(n)  c,
real(kind=8), dimension(n)  s 
)

Definition at line 5548 of file suews_util_minpack.f95.

Referenced by lmstr().

5548 
5549 !*****************************************************************************80
5550 !
5551 !! RWUPDT computes the decomposition of triangular matrix augmented by one row.
5552 !
5553 ! Discussion:
5554 !
5555 ! Given an N by N upper triangular matrix R, this subroutine
5556 ! computes the QR decomposition of the matrix formed when a row
5557 ! is added to R. If the row is specified by the vector W, then
5558 ! RWUPDT determines an orthogonal matrix Q such that when the
5559 ! N+1 by N matrix composed of R augmented by W is premultiplied
5560 ! by Q', the resulting matrix is upper trapezoidal.
5561 ! The matrix Q' is the product of N transformations
5562 !
5563 ! G(N)*G(N-1)* ... *G(1)
5564 !
5565 ! where G(I) is a Givens rotation in the (I,N+1) plane which eliminates
5566 ! elements in the (N+1)-st plane. RWUPDT also computes the product
5567 ! Q'*C where C is the (N+1)-vector (B,ALPHA). Q itself is not
5568 ! accumulated, rather the information to recover the G rotations is
5569 ! supplied.
5570 !
5571 ! Licensing:
5572 !
5573 ! This code is distributed under the GNU LGPL license.
5574 !
5575 ! Modified:
5576 !
5577 ! 06 April 2010
5578 !
5579 ! Author:
5580 !
5581 ! Original FORTRAN77 version by Jorge More, Burton Garbow, Kenneth Hillstrom.
5582 ! FORTRAN90 version by John Burkardt.
5583 !
5584 ! Reference:
5585 !
5586 ! Jorge More, Burton Garbow, Kenneth Hillstrom,
5587 ! User Guide for MINPACK-1,
5588 ! Technical Report ANL-80-74,
5589 ! Argonne National Laboratory, 1980.
5590 !
5591 ! Parameters:
5592 !
5593 ! Input, integer ( kind = 4 ) N, the order of R.
5594 !
5595 ! Input/output, real ( kind = 8 ) R(LDR,N). On input the upper triangular
5596 ! part of R must contain the matrix to be updated. On output R contains the
5597 ! updated triangular matrix.
5598 !
5599 ! Input, integer ( kind = 4 ) LDR, the leading dimension of the array R.
5600 ! LDR must not be less than N.
5601 !
5602 ! Input, real ( kind = 8 ) W(N), the row vector to be added to R.
5603 !
5604 ! Input/output, real ( kind = 8 ) B(N). On input, the first N elements
5605 ! of the vector C. On output the first N elements of the vector Q'*C.
5606 !
5607 ! Input/output, real ( kind = 8 ) ALPHA. On input, the (N+1)-st element
5608 ! of the vector C. On output the (N+1)-st element of the vector Q'*C.
5609 !
5610 ! Output, real ( kind = 8 ) C(N), S(N), the cosines and sines of the
5611 ! transforming Givens rotations.
5612 !
5613  implicit none
5614 
5615  integer(kind=4) ldr
5616  integer(kind=4) n
5617 
5618  real(kind=8) alpha
5619  real(kind=8) b(n)
5620  real(kind=8) c(n)
5621  real(kind=8) cotan
5622  integer(kind=4) i
5623  integer(kind=4) j
5624  real(kind=8) r(ldr, n)
5625  real(kind=8) rowj
5626  real(kind=8) s(n)
5627  real(kind=8) tan
5628  real(kind=8) temp
5629  real(kind=8) w(n)
5630 
5631  do j = 1, n
5632 
5633  rowj = w(j)
5634 !
5635 ! Apply the previous transformations to R(I,J), I=1,2,...,J-1, and to W(J).
5636 !
5637  do i = 1, j - 1
5638  temp = c(i)*r(i, j) + s(i)*rowj
5639  rowj = -s(i)*r(i, j) + c(i)*rowj
5640  r(i, j) = temp
5641  end do
5642 !
5643 ! Determine a Givens rotation which eliminates W(J).
5644 !
5645  c(j) = 1.0d+00
5646  s(j) = 0.0d+00
5647 
5648  if (rowj /= 0.0d+00) then
5649 
5650  if (abs(r(j, j)) < abs(rowj)) then
5651  cotan = r(j, j)/rowj
5652  s(j) = 0.5d+00/sqrt(0.25d+00 + 0.25d+00*cotan**2)
5653  c(j) = s(j)*cotan
5654  else
5655  tan = rowj/r(j, j)
5656  c(j) = 0.5d+00/sqrt(0.25d+00 + 0.25d+00*tan**2)
5657  s(j) = c(j)*tan
5658  end if
5659 !
5660 ! Apply the current transformation to R(J,J), B(J), and ALPHA.
5661 !
5662  r(j, j) = c(j)*r(j, j) + s(j)*rowj
5663  temp = c(j)*b(j) + s(j)*alpha
5664  alpha = -s(j)*b(j) + c(j)*alpha
5665  b(j) = temp
5666 
5667  end if
5668 
5669  end do
5670 
5671  return
Here is the caller graph for this function:

◆ timestamp()

subroutine timestamp ( )

Definition at line 5674 of file suews_util_minpack.f95.

5674 
5675 !*****************************************************************************80
5676 !
5677 !! TIMESTAMP prints the current YMDHMS date as a time stamp.
5678 !
5679 ! Example:
5680 !
5681 ! 31 May 2001 9:45:54.872 AM
5682 !
5683 ! Licensing:
5684 !
5685 ! This code is distributed under the GNU LGPL license.
5686 !
5687 ! Modified:
5688 !
5689 ! 18 May 2013
5690 !
5691 ! Author:
5692 !
5693 ! John Burkardt
5694 !
5695 ! Parameters:
5696 !
5697 ! None
5698 !
5699  implicit none
5700 
5701  character(len=8) ampm
5702  integer(kind=4) d
5703  integer(kind=4) h
5704  integer(kind=4) m
5705  integer(kind=4) mm
5706  character(len=9), parameter, dimension(12) :: month = (/ &
5707  'January ', 'February ', 'March ', 'April ', &
5708  'May ', 'June ', 'July ', 'August ', &
5709  'September', 'October ', 'November ', 'December '/)
5710  integer(kind=4) n
5711  integer(kind=4) s
5712  integer(kind=4) values(8)
5713  integer(kind=4) y
5714 
5715  call date_and_time(values=values)
5716 
5717  y = values(1)
5718  m = values(2)
5719  d = values(3)
5720  h = values(5)
5721  n = values(6)
5722  s = values(7)
5723  mm = values(8)
5724 
5725  if (h < 12) then
5726  ampm = 'AM'
5727  else if (h == 12) then
5728  if (n == 0 .and. s == 0) then
5729  ampm = 'Noon'
5730  else
5731  ampm = 'PM'
5732  end if
5733  else
5734  h = h - 12
5735  if (h < 12) then
5736  ampm = 'PM'
5737  else if (h == 12) then
5738  if (n == 0 .and. s == 0) then
5739  ampm = 'Midnight'
5740  else
5741  ampm = 'AM'
5742  end if
5743  end if
5744  end if
5745 
5746  write (*, '(i2.2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)') &
5747  d, trim(month(m)), y, h, ':', n, ':', s, '.', mm, trim(ampm)
5748 
5749  return
real(kind(1d0)) h