TABLE OF CONTENTS


ABINIT/m_levenberg_marquardt [ Modules ]

[ Top ] [ Modules ]

NAME

 m_levenberg_marquardt

FUNCTION

  Module containing routines for performing Levenberg-Marquardt
  nonlinear least-squares fit. These routines are the conversions
  to F90 of the Minpack F77 routines by Alan Miller

  TODO: For some reason, the routine lmder, which uses analytic
        evaluation of the Jacobian, did not work in testing (bug?).
        lmdif proved stable, but potentially uses more funtion
        evaluations.

COPYRIGHT

  Copyright (C) 2008-2018 ABINIT group (MS)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

28 #if defined HAVE_CONFIG_H
29 #include "config.h"
30 #endif
31 
32 #include "abi_common.h"
33 
34 module m_levenberg_marquardt
35 
36  use defs_basis
37  use m_abicore
38 ! MINPACK routines which are used by both LMDIF & LMDER
39 ! 25 October 2001:
40 !    Changed INTENT of iflag in several places to IN OUT.
41 !    Changed INTENT of fvec to IN OUT in user routine FCN.
42 !    Removed arguments diag and qtv from LMDIF & LMDER.
43 !    Replaced several DO loops with array operations.
44 ! amiller @ bigpond.net.au
45 
46  implicit none
47 
48  private
49 
50  public :: lmdif1
51  public :: lmdif
52  public :: lm_fit_print_info
53 
54 CONTAINS  !==============================================================================

m_levenberg_marquardt/enorm [ Functions ]

[ Top ] [ m_levenberg_marquardt ] [ Functions ]

NAME

FUNCTION

INPUTS

OUTPUT

NOTES

PARENTS

CHILDREN

SOURCE

1381  function enorm(n,x) result(fn_val)
1382 
1383  ! Code converted using TO_F90 by Alan Miller
1384  ! Date: 1999-12-09  Time: 12:45:34
1385 
1386 !Arguments ------------------------------------
1387 !scalars
1388 
1389 !This section has been created automatically by the script Abilint (TD).
1390 !Do not modify the following lines by hand.
1391 #undef ABI_FUNC
1392 #define ABI_FUNC 'enorm'
1393 !End of the abilint section
1394 
1395  real (dp)              :: fn_val
1396  integer, intent(in)    :: n
1397 !arrays
1398  real (dp), intent(in)  :: x(n)
1399 
1400 !Local variables-------------------------------
1401  integer   :: i
1402  real (dp) :: agiant, floatn, s1, s2, s3, xabs, x1max, x3max
1403  real (dp), parameter :: one = 1.0_dp, zero = 0.0_dp, rdwarf = 3.834e-20_dp, rgiant = 1.304e+19_dp
1404 ! *********************************************************************
1405 
1406  !  function enorm
1407 
1408  !  given an n-vector x, this function calculates the euclidean norm of x.
1409 
1410  !  the euclidean norm is computed by accumulating the sum of squares in
1411  !  three different sums.  The sums of squares for the small and large
1412  !  components are scaled so that no overflows occur.  Non-destructive
1413  !  underflows are permitted.  Underflows and overflows do not occur in the
1414  !  computation of the unscaled sum of squares for the intermediate
1415  !  components.  The definitions of small, intermediate and large components
1416  !  depend on two constants, rdwarf and rgiant.  The main restrictions on
1417  !  these constants are that rdwarf**2 not underflow and rgiant**2 not
1418  !  overflow.  The constants given here are suitable for every known computer.
1419 
1420  !  the function statement is
1421 
1422  !    REAL (dp) function enorm(n,x)
1423 
1424  !  where
1425 
1426  !    n is a positive integer input variable.
1427 
1428  !    x is an input array of length n.
1429 
1430  !  subprograms called
1431 
1432  !    fortran-supplied ... ABS,SQRT
1433 
1434  !  argonne national laboratory. minpack project. march 1980.
1435  !  burton s. garbow, kenneth e. hillstrom, jorge j. more
1436 
1437  !  **********
1438 
1439 
1440  s1 = zero
1441  s2 = zero
1442  s3 = zero
1443  x1max = zero
1444  x3max = zero
1445  floatn = n
1446  agiant = rgiant/floatn
1447  DO  i = 1, n
1448    xabs = ABS(x(i))
1449    IF (xabs > rdwarf .AND. xabs < agiant) GO TO 70
1450    IF (xabs <= rdwarf) GO TO 30
1451 
1452  !              sum for large components.
1453 
1454    IF (xabs <= x1max) GO TO 10
1455    s1 = one + s1*(x1max/xabs)**2
1456    x1max = xabs
1457    GO TO 20
1458 
1459    10 s1 = s1 + (xabs/x1max)**2
1460 
1461    20 GO TO 60
1462 
1463  !              sum for small components.
1464 
1465    30 IF (xabs <= x3max) GO TO 40
1466    s3 = one + s3*(x3max/xabs)**2
1467    x3max = xabs
1468    GO TO 60
1469 
1470    40 IF (xabs /= zero) s3 = s3 + (xabs/x3max)**2
1471 
1472    60 CYCLE
1473 
1474  !           sum for intermediate components.
1475 
1476    70 s2 = s2 + xabs**2
1477  END DO
1478 
1479  !     calculation of norm.
1480 
1481  IF (s1 == zero) GO TO 100
1482  fn_val = x1max*SQRT(s1 + (s2/x1max)/x1max)
1483  GO TO 120
1484 
1485  100 IF (s2 == zero) GO TO 110
1486  IF (s2 >= x3max) fn_val = SQRT(s2*(one + (x3max/s2)*(x3max*s3)))
1487  IF (s2 < x3max) fn_val = SQRT(x3max*((s2/x3max) + (x3max*s3)))
1488  GO TO 120
1489 
1490  110 fn_val = x3max*SQRT(s3)
1491 
1492  120 RETURN
1493 
1494  !     last card of function enorm.
1495 
1496  END FUNCTION enorm

m_levenberg_marquardt/fdjac2 [ Functions ]

[ Top ] [ m_levenberg_marquardt ] [ Functions ]

NAME

FUNCTION

INPUTS

OUTPUT

NOTES

PARENTS

CHILDREN

SOURCE

1518  SUBROUTINE fdjac2(fcn, m, n, x, fvec, fjac, iflag, epsfcn)
1519 
1520  ! Code converted using TO_F90 by Alan Miller
1521  ! Date: 1999-12-09  Time: 12:45:44
1522 
1523  ! N.B. Arguments LDFJAC & WA have been removed.
1524 
1525 !This section has been created automatically by the script Abilint (TD).
1526 !Do not modify the following lines by hand.
1527 #undef ABI_FUNC
1528 #define ABI_FUNC 'fdjac2'
1529 !End of the abilint section
1530 
1531  integer, intent(in)        :: m
1532  integer, intent(in)        :: n
1533  real (dp), intent(in out)  :: x(n)
1534  real (dp), intent(in)      :: fvec(m)
1535  real (dp), intent(out)     :: fjac(:,:)    ! fjac(ldfjac,n)
1536  integer, intent(in out)    :: iflag
1537  real (dp), intent(in)      :: epsfcn
1538 
1539  interface
1540    subroutine fcn(m, n, x, fvec, iflag, y, nfqre, nfqim, bsign, csign, multi_x_exp)
1541      implicit none
1542      integer, parameter         :: dp = KIND(1.0d0)
1543      integer, intent(in)        :: m, n
1544      real (dp), intent(in)      :: x(n)
1545      real (dp), intent(in out)  :: fvec(m)
1546      integer, intent(in out)    :: iflag
1547      integer, optional, intent(in) :: nfqre,nfqim,bsign,csign,multi_x_exp
1548      real (dp), optional, intent(in)      :: y(m)
1549    end subroutine fcn
1550  end interface
1551 
1552  !  **********
1553 
1554  !  subroutine fdjac2
1555 
1556  !  this subroutine computes a forward-difference approximation
1557  !  to the m by n jacobian matrix associated with a specified
1558  !  problem of m functions in n variables.
1559 
1560  !  the subroutine statement is
1561 
1562  !    subroutine fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa)
1563 
1564  !  where
1565 
1566  !    fcn is the name of the user-supplied subroutine which calculates the
1567  !      functions.  fcn must be declared in an external statement in the user
1568  !      calling program, and should be written as follows.
1569 
1570  !      subroutine fcn(m,n,x,fvec,iflag)
1571  !      integer m,n,iflag
1572  !      REAL (dp) x(n),fvec(m)
1573  !      ----------
1574  !      calculate the functions at x and
1575  !      return this vector in fvec.
1576  !      ----------
1577  !      return
1578  !      end
1579 
1580  !      the value of iflag should not be changed by fcn unless
1581  !      the user wants to terminate execution of fdjac2.
1582  !      in this case set iflag to a negative integer.
1583 
1584  !    m is a positive integer input variable set to the number of functions.
1585 
1586  !    n is a positive integer input variable set to the number of variables.
1587  !      n must not exceed m.
1588 
1589  !    x is an input array of length n.
1590 
1591  !    fvec is an input array of length m which must contain the
1592  !      functions evaluated at x.
1593 
1594  !    fjac is an output m by n array which contains the
1595  !      approximation to the jacobian matrix evaluated at x.
1596 
1597  !    ldfjac is a positive integer input variable not less than m
1598  !      which specifies the leading dimension of the array fjac.
1599 
1600  !    iflag is an integer variable which can be used to terminate
1601  !      the execution of fdjac2.  see description of fcn.
1602 
1603  !    epsfcn is an input variable used in determining a suitable step length
1604  !      for the forward-difference approximation.  This approximation assumes
1605  !      that the relative errors in the functions are of the order of epsfcn.
1606  !      If epsfcn is less than the machine precision, it is assumed that the
1607  !      relative errors in the functions are of the order of the machine
1608  !      precision.
1609 
1610  !    wa is a work array of length m.
1611 
1612  !  subprograms called
1613 
1614  !    user-supplied ...... fcn
1615 
1616  !    minpack-supplied ... dpmpar
1617 
1618  !    fortran-supplied ... ABS,MAX,SQRT
1619 
1620  !  argonne national laboratory. minpack project. march 1980.
1621  !  burton s. garbow, kenneth e. hillstrom, jorge j. more
1622 
1623  !  **********
1624  integer   :: j
1625  real (dp) :: eps, epsmch, h, temp, wa(m)
1626  real (dp), parameter :: zero = 0.0_dp
1627 
1628  !     epsmch is the machine precision.
1629 
1630  epsmch = EPSILON(zero)
1631 
1632  eps = SQRT(MAX(epsfcn, epsmch))
1633  DO  j = 1, n
1634    temp = x(j)
1635    h = eps*ABS(temp)
1636    IF (h == zero) h = eps
1637    x(j) = temp + h
1638    CALL fcn(m, n, x, wa, iflag)
1639    IF (iflag < 0) EXIT
1640    x(j) = temp
1641    fjac(1:m,j) = (wa(1:m) - fvec(1:m))/h
1642  END DO
1643 
1644  RETURN
1645 
1646 !     last card of subroutine fdjac2.
1647 
1648  end subroutine fdjac2

m_levenberg_marquardt/lm_fit_print_info [ Functions ]

[ Top ] [ m_levenberg_marquardt ] [ Functions ]

NAME

FUNCTION

INPUTS

OUTPUT

NOTES

PARENTS

CHILDREN

SOURCE

1667 subroutine lm_fit_print_info(info,msg)
1668 
1669 
1670 !This section has been created automatically by the script Abilint (TD).
1671 !Do not modify the following lines by hand.
1672 #undef ABI_FUNC
1673 #define ABI_FUNC 'lm_fit_print_info'
1674 !End of the abilint section
1675 
1676   implicit none
1677 
1678   integer, intent(in) :: info
1679   character(len=500)  :: msg
1680 
1681   select case (info)
1682   case (:-1)
1683      write(msg,'(a,i0)') 'STOP LM fit: fcn in LM fit returned info = ', -info
1684   case (0)
1685      write(msg,'(a)') 'ERROR LM fit: improper values for input parameters in LM fit'
1686   case (1:3)
1687      write(msg,'(a)') 'LM fit converged'
1688   case (4)
1689      write(msg,'(a)') 'ERROR LM fit: residuals orthogonal to the jacobian, there may be an error in fcn.'
1690   case (5)
1691      write(msg,'(a)') 'ERROR LM fit: too many calls to fcn, either slow convergence, or an error in opt.'
1692   case (6:7)
1693      write(msg,'(a)') 'ERROR LM fit: tol was set too small'
1694   case default
1695      write(msg,'(a,i0)') 'ERROR LM fit: unknown info = ',info
1696   end select
1697 
1698   return
1699 
1700 end subroutine lm_fit_print_info
1701 
1702 end module m_levenberg_marquardt

m_levenberg_marquardt/lmdif [ Functions ]

[ Top ] [ m_levenberg_marquardt ] [ Functions ]

NAME

FUNCTION

INPUTS

OUTPUT

NOTES

PARENTS

CHILDREN

SOURCE

251  SUBROUTINE lmdif(fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn,  &
252                   mode, factor, nprint, info, nfev, fjac, ipvt)
253 
254  ! N.B. Arguments LDFJAC, DIAG, QTF, WA1, WA2, WA3 & WA4 have been removed.
255 
256 !This section has been created automatically by the script Abilint (TD).
257 !Do not modify the following lines by hand.
258 #undef ABI_FUNC
259 #define ABI_FUNC 'lmdif'
260 !End of the abilint section
261 
262  integer, intent(in)        :: m
263  integer, intent(in)        :: n
264  real (dp), intent(in out)  :: x(:)
265  real (dp), intent(out)     :: fvec(:)
266  real (dp), intent(in)      :: ftol
267  real (dp), intent(in)      :: xtol
268  real (dp), intent(in out)  :: gtol
269  integer, intent(in out)    :: maxfev
270  real (dp), intent(in out)  :: epsfcn
271  integer, intent(in)        :: mode
272  real (dp), intent(in)      :: factor
273  integer, intent(in)        :: nprint
274  integer, intent(out)       :: info
275  integer, intent(out)       :: nfev
276  real (dp), intent(out)     :: fjac(:,:)    ! fjac(ldfjac,n)
277  integer, intent(out)       :: ipvt(:)
278 
279  ! external fcn
280 
281  interface
282    subroutine fcn(m, n, x, fvec, iflag, y, nfqre, nfqim, bsign, csign, multi_x_exp)
283      implicit none
284      integer, parameter         :: dp = KIND(1.0d0)
285      integer, intent(in)        :: m, n
286      real (dp), intent(in)      :: x(n)
287      real (dp), intent(in out)  :: fvec(m)
288      integer, intent(in out)    :: iflag
289      integer, optional, intent(in) :: nfqre,nfqim,bsign,csign,multi_x_exp
290      real (dp), optional, intent(in)      :: y(m)
291    end subroutine fcn
292  end interface
293 
294  !  **********
295 
296  !  subroutine lmdif
297 
298  !  The purpose of lmdif is to minimize the sum of the squares of m nonlinear
299  !  functions in n variables by a modification of the Levenberg-Marquardt
300  !  algorithm.  The user must provide a subroutine which calculates the
301  !  functions.  The jacobian is then calculated by a forward-difference
302  !  approximation.
303 
304  !  the subroutine statement is
305 
306  !    subroutine lmdif(fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn,
307  !                     diag, mode, factor, nprint, info, nfev, fjac,
308  !                     ldfjac, ipvt, qtf, wa1, wa2, wa3, wa4)
309 
310  ! N.B. 7 of these arguments have been removed in this version.
311 
312  !  where
313 
314  !    fcn is the name of the user-supplied subroutine which calculates the
315  !      functions.  fcn must be declared in an external statement in the user
316  !      calling program, and should be written as follows.
317 
318  !      subroutine fcn(m, n, x, fvec, iflag)
319  !      integer m, n, iflag
320  !      REAL (dp) x(:), fvec(m)
321  !      ----------
322  !      calculate the functions at x and return this vector in fvec.
323  !      ----------
324  !      return
325  !      end
326 
327  !      the value of iflag should not be changed by fcn unless
328  !      the user wants to terminate execution of lmdif.
329  !      in this case set iflag to a negative integer.
330 
331  !    m is a positive integer input variable set to the number of functions.
332 
333  !    n is a positive integer input variable set to the number of variables.
334  !      n must not exceed m.
335 
336  !    x is an array of length n.  On input x must contain an initial estimate
337  !      of the solution vector.  On output x contains the final estimate of the
338  !      solution vector.
339 
340  !    fvec is an output array of length m which contains
341  !      the functions evaluated at the output x.
342 
343  !    ftol is a nonnegative input variable.  Termination occurs when both the
344  !      actual and predicted relative reductions in the sum of squares are at
345  !      most ftol.  Therefore, ftol measures the relative error desired
346  !      in the sum of squares.
347 
348  !    xtol is a nonnegative input variable.  Termination occurs when the
349  !      relative error between two consecutive iterates is at most xtol.
350  !      Therefore, xtol measures the relative error desired in the approximate
351  !      solution.
352 
353  !    gtol is a nonnegative input variable.  Termination occurs when the cosine
354  !      of the angle between fvec and any column of the jacobian is at most
355  !      gtol in absolute value.  Therefore, gtol measures the orthogonality
356  !      desired between the function vector and the columns of the jacobian.
357 
358  !    maxfev is a positive integer input variable.  Termination occurs when the
359  !      number of calls to fcn is at least maxfev by the end of an iteration.
360 
361  !    epsfcn is an input variable used in determining a suitable step length
362  !      for the forward-difference approximation.  This approximation assumes
363  !      that the relative errors in the functions are of the order of epsfcn.
364  !      If epsfcn is less than the machine precision, it is assumed that the
365  !      relative errors in the functions are of the order of the machine
366  !      precision.
367 
368  !    diag is an array of length n.  If mode = 1 (see below), diag is
369  !      internally set.  If mode = 2, diag must contain positive entries that
370  !      serve as multiplicative scale factors for the variables.
371 
372  !    mode is an integer input variable.  If mode = 1, the variables will be
373  !      scaled internally.  If mode = 2, the scaling is specified by the input
374  !      diag. other values of mode are equivalent to mode = 1.
375 
376  !    factor is a positive input variable used in determining the initial step
377  !      bound.  This bound is set to the product of factor and the euclidean
378  !      norm of diag*x if nonzero, or else to factor itself.  In most cases
379  !      factor should lie in the interval (.1,100.). 100. is a generally
380  !      recommended value.
381 
382  !    nprint is an integer input variable that enables controlled printing of
383  !      iterates if it is positive.  In this case, fcn is called with iflag = 0
384  !      at the beginning of the first iteration and every nprint iterations
385  !      thereafter and immediately prior to return, with x and fvec available
386  !      for printing.  If nprint is not positive, no special calls
387  !      of fcn with iflag = 0 are made.
388 
389  !    info is an integer output variable.  If the user has terminated
390  !      execution, info is set to the (negative) value of iflag.
391  !      See description of fcn.  Otherwise, info is set as follows.
392 
393  !      info = 0  improper input parameters.
394 
395  !      info = 1  both actual and predicted relative reductions
396  !                in the sum of squares are at most ftol.
397 
398  !      info = 2  relative error between two consecutive iterates <= xtol.
399 
400  !      info = 3  conditions for info = 1 and info = 2 both hold.
401 
402  !      info = 4  the cosine of the angle between fvec and any column of
403  !                the Jacobian is at most gtol in absolute value.
404 
405  !      info = 5  number of calls to fcn has reached or exceeded maxfev.
406 
407  !      info = 6  ftol is too small. no further reduction in
408  !                the sum of squares is possible.
409 
410  !      info = 7  xtol is too small. no further improvement in
411  !                the approximate solution x is possible.
412 
413  !      info = 8  gtol is too small. fvec is orthogonal to the
414  !                columns of the jacobian to machine precision.
415 
416  !    nfev is an integer output variable set to the number of calls to fcn.
417 
418  !    fjac is an output m by n array. the upper n by n submatrix
419  !      of fjac contains an upper triangular matrix r with
420  !      diagonal elements of nonincreasing magnitude such that
421 
422  !             t     t           t
423  !            p *(jac *jac)*p = r *r,
424 
425  !      where p is a permutation matrix and jac is the final calculated
426  !      Jacobian.  Column j of p is column ipvt(j) (see below) of the
427  !      identity matrix. the lower trapezoidal part of fjac contains
428  !      information generated during the computation of r.
429 
430  !    ldfjac is a positive integer input variable not less than m
431  !      which specifies the leading dimension of the array fjac.
432 
433  !    ipvt is an integer output array of length n.  ipvt defines a permutation
434  !      matrix p such that jac*p = q*r, where jac is the final calculated
435  !      jacobian, q is orthogonal (not stored), and r is upper triangular
436  !      with diagonal elements of nonincreasing magnitude.
437  !      Column j of p is column ipvt(j) of the identity matrix.
438 
439  !    qtf is an output array of length n which contains
440  !      the first n elements of the vector (q transpose)*fvec.
441 
442  !    wa1, wa2, and wa3 are work arrays of length n.
443 
444  !    wa4 is a work array of length m.
445 
446  !  subprograms called
447 
448  !    user-supplied ...... fcn
449 
450  !    minpack-supplied ... dpmpar,enorm,fdjac2,lmpar,qrfac
451 
452  !    fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod
453 
454  !  argonne national laboratory. minpack project. march 1980.
455  !  burton s. garbow, kenneth e. hillstrom, jorge j. more
456 
457  !  **********
458  integer   :: i, iflag, iter, j, l
459  real (dp) :: actred, delta, dirder, epsmch, fnorm, fnorm1, gnorm
460  real(dp)  :: par, pnorm, prered, ratio, sum, temp, temp1, temp2, xnorm
461  real (dp) :: diag(n), qtf(n), wa1(n), wa2(n), wa3(n), wa4(m)
462  real (dp), parameter :: one = 1.0_dp, p1 = 0.1_dp, p5 = 0.5_dp
463  real(dp),parameter ::   p25 = 0.25_dp, p75 = 0.75_dp, p0001 = 0.0001_dp, zero = 0.0_dp
464 
465  !     epsmch is the machine precision.
466 
467  epsmch = EPSILON(zero)
468 
469  info = 0
470  iflag = 0
471  nfev = 0
472 
473  !     check the input parameters for errors.
474 
475  IF (n <= 0 .OR. m < n .OR. ftol < zero .OR. xtol < zero .OR. gtol < zero  &
476 &     .OR. maxfev <= 0 .OR. factor <= zero) GO TO 300
477  IF (mode /= 2) GO TO 20
478  DO  j = 1, n
479    IF (diag(j) <= zero) GO TO 300
480  END DO
481 
482  !     evaluate the function at the starting point and calculate its norm.
483 
484  20 iflag = 1
485  CALL fcn(m, n, x, fvec, iflag)
486  nfev = 1
487  IF (iflag < 0) GO TO 300
488  fnorm = enorm(m, fvec)
489 
490  !     initialize levenberg-marquardt parameter and iteration counter.
491 
492  par = zero
493  iter = 1
494 
495  !     beginning of the outer loop.
496 
497  !        calculate the jacobian matrix.
498 
499  30 iflag = 2
500  CALL fdjac2(fcn, m, n, x, fvec, fjac, iflag, epsfcn)
501  nfev = nfev + n
502  IF (iflag < 0) GO TO 300
503 
504  !        If requested, call fcn to enable printing of iterates.
505 
506  IF (nprint <= 0) GO TO 40
507  iflag = 0
508  IF (MOD(iter-1,nprint) == 0) THEN
509    CALL fcn(m, n, x, fvec, iflag)
510  END IF
511  IF (iflag < 0) GO TO 300
512 
513  !        Compute the qr factorization of the jacobian.
514 
515  40 CALL qrfac(m, n, fjac, .true., ipvt, wa1, wa2)
516 
517  !        On the first iteration and if mode is 1, scale according
518  !        to the norms of the columns of the initial jacobian.
519 
520  IF (iter /= 1) GO TO 80
521  IF (mode == 2) GO TO 60
522  DO  j = 1, n
523    diag(j) = wa2(j)
524    IF (wa2(j) == zero) diag(j) = one
525  END DO
526 
527  !        On the first iteration, calculate the norm of the scaled x
528  !        and initialize the step bound delta.
529 
530  60 wa3(1:n) = diag(1:n)*x(1:n)
531  xnorm = enorm(n, wa3)
532  delta = factor*xnorm
533  IF (delta == zero) delta = factor
534 
535  !        Form (q transpose)*fvec and store the first n components in qtf.
536 
537  80 wa4(1:m) = fvec(1:m)
538  DO  j = 1, n
539    IF (fjac(j,j) == zero) GO TO 120
540    sum = DOT_PRODUCT( fjac(j:m,j), wa4(j:m) )
541    temp = -sum/fjac(j,j)
542    DO  i = j, m
543      wa4(i) = wa4(i) + fjac(i,j)*temp
544    END DO
545    120 fjac(j,j) = wa1(j)
546    qtf(j) = wa4(j)
547  END DO
548 
549  !        compute the norm of the scaled gradient.
550 
551  gnorm = zero
552  IF (fnorm == zero) GO TO 170
553  DO  j = 1, n
554    l = ipvt(j)
555    IF (wa2(l) == zero) CYCLE
556    sum = zero
557    DO  i = 1, j
558      sum = sum + fjac(i,j)*(qtf(i)/fnorm)
559    END DO
560    gnorm = MAX(gnorm, ABS(sum/wa2(l)))
561  END DO
562 
563  !        test for convergence of the gradient norm.
564 
565  170 IF (gnorm <= gtol) info = 4
566  IF (info /= 0) GO TO 300
567 
568  !        rescale if necessary.
569 
570  IF (mode == 2) GO TO 200
571  DO  j = 1, n
572    diag(j) = MAX(diag(j), wa2(j))
573  END DO
574 
575  !        beginning of the inner loop.
576 
577  !           determine the Levenberg-Marquardt parameter.
578 
579  200 CALL lmpar(n, fjac, ipvt, diag, qtf, delta, par, wa1, wa2)
580 
581  !           store the direction p and x + p. calculate the norm of p.
582 
583  DO  j = 1, n
584    wa1(j) = -wa1(j)
585    wa2(j) = x(j) + wa1(j)
586    wa3(j) = diag(j)*wa1(j)
587  END DO
588  pnorm = enorm(n, wa3)
589 
590  !           on the first iteration, adjust the initial step bound.
591 
592  IF (iter == 1) delta = MIN(delta, pnorm)
593 
594  !           evaluate the function at x + p and calculate its norm.
595 
596  iflag = 1
597  CALL fcn(m, n, wa2, wa4, iflag)
598  nfev = nfev + 1
599  IF (iflag < 0) GO TO 300
600  fnorm1 = enorm(m, wa4)
601 
602  !           compute the scaled actual reduction.
603 
604  actred = -one
605  IF (p1*fnorm1 < fnorm) actred = one - (fnorm1/fnorm)**2
606 
607  !           Compute the scaled predicted reduction and
608  !           the scaled directional derivative.
609 
610  DO  j = 1, n
611    wa3(j) = zero
612    l = ipvt(j)
613    temp = wa1(l)
614    DO  i = 1, j
615      wa3(i) = wa3(i) + fjac(i,j)*temp
616    END DO
617  END DO
618  temp1 = enorm(n,wa3)/fnorm
619  temp2 = (SQRT(par)*pnorm)/fnorm
620  prered = temp1**2 + temp2**2/p5
621  dirder = -(temp1**2 + temp2**2)
622 
623  !           compute the ratio of the actual to the predicted reduction.
624 
625  ratio = zero
626  IF (prered /= zero) ratio = actred/prered
627 
628  !           update the step bound.
629 
630  IF (ratio <= p25) THEN
631    IF (actred >= zero) temp = p5
632    IF (actred < zero) temp = p5*dirder/(dirder + p5*actred)
633    IF (p1*fnorm1 >= fnorm .OR. temp < p1) temp = p1
634    delta = temp*MIN(delta,pnorm/p1)
635    par = par/temp
636  ELSE
637    IF (par /= zero .AND. ratio < p75) GO TO 260
638    delta = pnorm/p5
639    par = p5*par
640  END IF
641 
642  !           test for successful iteration.
643 
644  260 IF (ratio < p0001) GO TO 290
645 
646  !           successful iteration. update x, fvec, and their norms.
647 
648  DO  j = 1, n
649    x(j) = wa2(j)
650    wa2(j) = diag(j)*x(j)
651  END DO
652  fvec(1:m) = wa4(1:m)
653  xnorm = enorm(n, wa2)
654  fnorm = fnorm1
655  iter = iter + 1
656 
657  !           tests for convergence.
658 
659  290 IF (ABS(actred) <= ftol .AND. prered <= ftol .AND. p5*ratio <= one) info = 1
660  IF (delta <= xtol*xnorm) info = 2
661  IF (ABS(actred) <= ftol .AND. prered <= ftol  &
662 &     .AND. p5*ratio <= one .AND. info == 2) info = 3
663  IF (info /= 0) GO TO 300
664 
665  !           tests for termination and stringent tolerances.
666 
667  IF (nfev >= maxfev) info = 5
668  IF (ABS(actred) <= epsmch .AND. prered <= epsmch  &
669 &     .AND. p5*ratio <= one) info = 6
670  IF (delta <= epsmch*xnorm) info = 7
671  IF (gnorm <= epsmch) info = 8
672  IF (info /= 0) GO TO 300
673 
674  !           end of the inner loop. repeat if iteration unsuccessful.
675 
676  IF (ratio < p0001) GO TO 200
677 
678  !        end of the outer loop.
679 
680  GO TO 30
681 
682  !     termination, either normal or user imposed.
683 
684  300 IF (iflag < 0) info = iflag
685  iflag = 0
686  IF (nprint > 0) THEN
687    CALL fcn(m, n, x, fvec, iflag)
688  END IF
689  RETURN
690 
691  !     last card of subroutine lmdif.
692 
693  end subroutine lmdif

m_levenberg_marquardt/lmdif1 [ Functions ]

[ Top ] [ m_levenberg_marquardt ] [ Functions ]

NAME

FUNCTION

INPUTS

OUTPUT

NOTES

PARENTS

CHILDREN

SOURCE

 74  subroutine lmdif1(fcn, m, n, x, fvec, tol, info, iwa)
 75 
 76  ! code converted using to_f90 by alan miller
 77  ! date: 1999-12-11  time: 00:51:44
 78 
 79  ! n.b. arguments wa & lwa have been removed.
 80 
 81 !This section has been created automatically by the script Abilint (TD).
 82 !Do not modify the following lines by hand.
 83 #undef ABI_FUNC
 84 #define ABI_FUNC 'lmdif1'
 85 !End of the abilint section
 86 
 87  integer, intent(in)        :: m
 88  integer, intent(in)        :: n
 89  real (dp), intent(in out)  :: x(:)
 90  real (dp), intent(out)     :: fvec(:)
 91  real (dp), intent(in)      :: tol
 92  integer, intent(out)       :: info
 93  integer, intent(out)       :: iwa(:)
 94 
 95  ! external fcn
 96 
 97  interface
 98    subroutine fcn(m, n, x, fvec, iflag, y, nfqre, nfqim, bsign, csign, multi_x_exp)
 99      implicit none
100      integer, parameter         :: dp = KIND(1.0d0)
101      integer, intent(in)        :: m, n
102      real (dp), intent(in)      :: x(n)
103      real (dp), intent(in out)  :: fvec(m)
104      integer, intent(in out)    :: iflag
105      integer, optional, intent(in) :: nfqre,nfqim,bsign,csign,multi_x_exp
106      real (dp), optional, intent(in)      :: y(m)
107    end subroutine fcn
108  end interface
109 
110  !  **********
111 
112  !  subroutine lmdif1
113 
114  !  The purpose of lmdif1 is to minimize the sum of the squares of m nonlinear
115  !  functions in n variables by a modification of the Levenberg-Marquardt
116  !  algorithm.  This is done by using the more general least-squares
117  !  solver lmdif.  The user must provide a subroutine which calculates the
118  !  functions.  The jacobian is then calculated by a forward-difference
119  !  approximation.
120 
121  !  the subroutine statement is
122 
123  !    subroutine lmdif1(fcn, m, n, x, fvec, tol, info, iwa)
124 
125  !  where
126 
127  !    fcn is the name of the user-supplied subroutine which calculates
128  !      the functions.  fcn must be declared in an external statement in the
129  !      user calling program, and should be written as follows.
130 
131  !      subroutine fcn(m, n, x, fvec, iflag)
132  !      integer m, n, iflag
133  !      REAL (dp) x(n), fvec(m)
134  !      ----------
135  !      calculate the functions at x and return this vector in fvec.
136  !      ----------
137  !      return
138  !      end
139 
140  !      the value of iflag should not be changed by fcn unless
141  !      the user wants to terminate execution of lmdif1.
142  !      In this case set iflag to a negative integer.
143 
144  !    m is a positive integer input variable set to the number of functions.
145 
146  !    n is a positive integer input variable set to the number of variables.
147  !      n must not exceed m.
148 
149  !    x is an array of length n.  On input x must contain an initial estimate
150  !      of the solution vector.  On output x contains the final estimate of
151  !      the solution vector.
152 
153  !    fvec is an output array of length m which contains
154  !      the functions evaluated at the output x.
155 
156  !    tol is a nonnegative input variable.  Termination occurs when the
157  !      algorithm estimates either that the relative error in the sum of
158  !      squares is at most tol or that the relative error between x and the
159  !      solution is at most tol.
160 
161  !    info is an integer output variable.  If the user has terminated execution,
162  !      info is set to the (negative) value of iflag.  See description of fcn.
163  !      Otherwise, info is set as follows.
164 
165  !      info = 0  improper input parameters.
166 
167  !      info = 1  algorithm estimates that the relative error
168  !                in the sum of squares is at most tol.
169 
170  !      info = 2  algorithm estimates that the relative error
171  !                between x and the solution is at most tol.
172 
173  !      info = 3  conditions for info = 1 and info = 2 both hold.
174 
175  !      info = 4  fvec is orthogonal to the columns of the
176  !                jacobian to machine precision.
177 
178  !      info = 5  number of calls to fcn has reached or exceeded 200*(n+1).
179 
180  !      info = 6  tol is too small. no further reduction in
181  !                the sum of squares is possible.
182 
183  !      info = 7  tol is too small.  No further improvement in
184  !                the approximate solution x is possible.
185 
186  !    iwa is an integer work array of length n.
187 
188  !    wa is a work array of length lwa.
189 
190  !    lwa is a positive integer input variable not less than m*n+5*n+m.
191 
192  !  subprograms called
193 
194  !    user-supplied ...... fcn
195 
196  !    minpack-supplied ... lmdif
197 
198  !  argonne national laboratory. minpack project. march 1980.
199  !  burton s. garbow, kenneth e. hillstrom, jorge j. more
200 
201  !  **********
202  integer   :: maxfev, mode, nfev, nprint
203  real (dp) :: epsfcn, ftol, gtol, xtol, fjac(m,n)
204  real (dp), parameter :: factor = 100._dp , zero = 0.0_dp
205 
206  info = 0
207 
208  !     check the input parameters for errors.
209 
210  IF (n <= 0 .OR. m < n .OR. tol < zero) GO TO 10
211 
212  !     call lmdif.
213 
214  maxfev = 200*(n + 1)
215  ftol = tol
216  xtol = tol
217  gtol = zero
218  epsfcn = zero
219  mode = 1
220  nprint = 0
221  CALL lmdif(fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn,   &
222 &           mode, factor, nprint, info, nfev, fjac, iwa)
223  IF (info == 8) info = 4
224 
225  10 RETURN
226 
227  !     last card of subroutine lmdif1.
228 
229  end subroutine lmdif1

m_levenberg_marquardt/lmpar [ Functions ]

[ Top ] [ m_levenberg_marquardt ] [ Functions ]

NAME

FUNCTION

INPUTS

OUTPUT

NOTES

PARENTS

CHILDREN

SOURCE

715  subroutine lmpar(n, r, ipvt, diag, qtb, delta, par, x, sdiag)
716 
717  ! Code converted using TO_F90 by Alan Miller
718  ! Date: 1999-12-09  Time: 12:46:12
719 
720  ! N.B. Arguments LDR, WA1 & WA2 have been removed.
721 
722 !This section has been created automatically by the script Abilint (TD).
723 !Do not modify the following lines by hand.
724 #undef ABI_FUNC
725 #define ABI_FUNC 'lmpar'
726 !End of the abilint section
727 
728  integer, intent(in)        :: n
729  real (dp), intent(in out)  :: r(:,:)
730  integer, intent(in)        :: ipvt(:)
731  real (dp), intent(in)      :: diag(:)
732  real (dp), intent(in)      :: qtb(:)
733  real (dp), intent(in)      :: delta
734  real (dp), intent(out)     :: par
735  real (dp), intent(out)     :: x(:)
736  real (dp), intent(out)     :: sdiag(:)
737 
738  !  **********
739 
740  !  subroutine lmpar
741 
742  !  Given an m by n matrix a, an n by n nonsingular diagonal matrix d,
743  !  an m-vector b, and a positive number delta, the problem is to determine a
744  !  value for the parameter par such that if x solves the system
745 
746  !        a*x = b ,     sqrt(par)*d*x = 0 ,
747 
748  !  in the least squares sense, and dxnorm is the euclidean
749  !  norm of d*x, then either par is zero and
750 
751  !        (dxnorm-delta) <= 0.1*delta ,
752 
753  !  or par is positive and
754 
755  !        abs(dxnorm-delta) <= 0.1*delta .
756 
757  !  This subroutine completes the solution of the problem if it is provided
758  !  with the necessary information from the r factorization, with column
759  !  qpivoting, of a.  That is, if a*p = q*r, where p is a permutation matrix,
760  !  q has orthogonal columns, and r is an upper triangular matrix with diagonal
761  !  elements of nonincreasing magnitude, then lmpar expects the full upper
762  !  triangle of r, the permutation matrix p, and the first n components of
763  !  (q transpose)*b.
764  !  On output lmpar also provides an upper triangular matrix s such that
765 
766  !         t   t                   t
767  !        p *(a *a + par*d*d)*p = s *s .
768 
769  !  s is employed within lmpar and may be of separate interest.
770 
771  !  Only a few iterations are generally needed for convergence of the algorithm.
772  !  If, however, the limit of 10 iterations is reached, then the output par
773  !  will contain the best value obtained so far.
774 
775  !  the subroutine statement is
776 
777  !    subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag, wa1,wa2)
778 
779  !  where
780 
781  !    n is a positive integer input variable set to the order of r.
782 
783  !    r is an n by n array. on input the full upper triangle
784  !      must contain the full upper triangle of the matrix r.
785  !      On output the full upper triangle is unaltered, and the
786  !      strict lower triangle contains the strict upper triangle
787  !      (transposed) of the upper triangular matrix s.
788 
789  !    ldr is a positive integer input variable not less than n
790  !      which specifies the leading dimension of the array r.
791 
792  !    ipvt is an integer input array of length n which defines the
793  !      permutation matrix p such that a*p = q*r. column j of p
794  !      is column ipvt(j) of the identity matrix.
795 
796  !    diag is an input array of length n which must contain the
797  !      diagonal elements of the matrix d.
798 
799  !    qtb is an input array of length n which must contain the first
800  !      n elements of the vector (q transpose)*b.
801 
802  !    delta is a positive input variable which specifies an upper
803  !      bound on the euclidean norm of d*x.
804 
805  !    par is a nonnegative variable. on input par contains an
806  !      initial estimate of the levenberg-marquardt parameter.
807  !      on output par contains the final estimate.
808 
809  !    x is an output array of length n which contains the least
810  !      squares solution of the system a*x = b, sqrt(par)*d*x = 0,
811  !      for the output par.
812 
813  !    sdiag is an output array of length n which contains the
814  !      diagonal elements of the upper triangular matrix s.
815 
816  !    wa1 and wa2 are work arrays of length n.
817 
818  !  subprograms called
819 
820  !    minpack-supplied ... dpmpar,enorm,qrsolv
821 
822  !    fortran-supplied ... ABS,MAX,MIN,SQRT
823 
824  !  argonne national laboratory. minpack project. march 1980.
825  !  burton s. garbow, kenneth e. hillstrom, jorge j. more
826 
827  !  **********
828  integer   :: iter, j, jm1, jp1, k, l, nsing
829  real (dp) :: dxnorm, dwarf, fp, gnorm, parc, parl, paru, sum, temp
830  real (dp) :: wa1(n), wa2(n)
831  real (dp), parameter :: p1 = 0.1_dp, p001 = 0.001_dp, zero = 0.0_dp
832 
833  !     dwarf is the smallest positive magnitude.
834 
835  dwarf = TINY(zero)
836 
837  !     compute and store in x the gauss-newton direction. if the
838  !     jacobian is rank-deficient, obtain a least squares solution.
839 
840  nsing = n
841  DO  j = 1, n
842    wa1(j) = qtb(j)
843    IF (r(j,j) == zero .AND. nsing == n) nsing = j - 1
844    IF (nsing < n) wa1(j) = zero
845  END DO
846 
847  DO  k = 1, nsing
848    j = nsing - k + 1
849    wa1(j) = wa1(j)/r(j,j)
850    temp = wa1(j)
851    jm1 = j - 1
852    wa1(1:jm1) = wa1(1:jm1) - r(1:jm1,j)*temp
853  END DO
854 
855  DO  j = 1, n
856    l = ipvt(j)
857    x(l) = wa1(j)
858  END DO
859 
860  !     initialize the iteration counter.
861  !     evaluate the function at the origin, and test
862  !     for acceptance of the gauss-newton direction.
863 
864  iter = 0
865  wa2(1:n) = diag(1:n)*x(1:n)
866  dxnorm = enorm(n, wa2)
867  fp = dxnorm - delta
868  IF (fp <= p1*delta) GO TO 220
869 
870  !     if the jacobian is not rank deficient, the newton
871  !     step provides a lower bound, parl, for the zero of
872  !     the function.  Otherwise set this bound to zero.
873 
874  parl = zero
875  IF (nsing < n) GO TO 120
876  DO  j = 1, n
877    l = ipvt(j)
878    wa1(j) = diag(l)*(wa2(l)/dxnorm)
879  END DO
880  DO  j = 1, n
881    sum = DOT_PRODUCT( r(1:j-1,j), wa1(1:j-1) )
882    wa1(j) = (wa1(j) - sum)/r(j,j)
883  END DO
884  temp = enorm(n,wa1)
885  parl = ((fp/delta)/temp)/temp
886 
887  !     calculate an upper bound, paru, for the zero of the function.
888 
889  120 DO  j = 1, n
890    sum = DOT_PRODUCT( r(1:j,j), qtb(1:j) )
891    l = ipvt(j)
892    wa1(j) = sum/diag(l)
893  END DO
894  gnorm = enorm(n,wa1)
895  paru = gnorm/delta
896  IF (paru == zero) paru = dwarf/MIN(delta,p1)
897 
898  !     if the input par lies outside of the interval (parl,paru),
899  !     set par to the closer endpoint.
900 
901  par = MAX(par,parl)
902  par = MIN(par,paru)
903  IF (par == zero) par = gnorm/dxnorm
904 
905  !     beginning of an iteration.
906 
907  150 iter = iter + 1
908 
909  !        evaluate the function at the current value of par.
910 
911  IF (par == zero) par = MAX(dwarf, p001*paru)
912  temp = SQRT(par)
913  wa1(1:n) = temp*diag(1:n)
914  CALL qrsolv(n, r, ipvt, wa1, qtb, x, sdiag)
915  wa2(1:n) = diag(1:n)*x(1:n)
916  dxnorm = enorm(n, wa2)
917  temp = fp
918  fp = dxnorm - delta
919 
920  !        if the function is small enough, accept the current value
921  !        of par. also test for the exceptional cases where parl
922  !        is zero or the number of iterations has reached 10.
923 
924  IF (ABS(fp) <= p1*delta .OR. parl == zero .AND. fp <= temp  &
925 &     .AND. temp < zero .OR. iter == 10) GO TO 220
926 
927  !        compute the newton correction.
928 
929  DO  j = 1, n
930    l = ipvt(j)
931    wa1(j) = diag(l)*(wa2(l)/dxnorm)
932  END DO
933  DO  j = 1, n
934    wa1(j) = wa1(j)/sdiag(j)
935    temp = wa1(j)
936    jp1 = j + 1
937    wa1(jp1:n) = wa1(jp1:n) - r(jp1:n,j)*temp
938  END DO
939  temp = enorm(n,wa1)
940  parc = ((fp/delta)/temp)/temp
941 
942  !        depending on the sign of the function, update parl or paru.
943 
944  IF (fp > zero) parl = MAX(parl,par)
945  IF (fp < zero) paru = MIN(paru,par)
946 
947  !        compute an improved estimate for par.
948 
949  par = MAX(parl, par+parc)
950 
951  !        end of an iteration.
952 
953  GO TO 150
954 
955  !     termination.
956 
957  220 IF (iter == 0) par = zero
958  RETURN
959 
960  !     last card of subroutine lmpar.
961 
962  end subroutine lmpar

m_levenberg_marquardt/qrfac [ Functions ]

[ Top ] [ m_levenberg_marquardt ] [ Functions ]

NAME

FUNCTION

INPUTS

OUTPUT

NOTES

PARENTS

CHILDREN

SOURCE

 984  subroutine qrfac(m, n, a, pivot, ipvt, rdiag, acnorm)
 985 
 986  ! Code converted using TO_F90 by Alan Miller
 987  ! Date: 1999-12-09  Time: 12:46:17
 988 
 989  ! N.B. Arguments LDA, LIPVT & WA have been removed.
 990 
 991 !This section has been created automatically by the script Abilint (TD).
 992 !Do not modify the following lines by hand.
 993 #undef ABI_FUNC
 994 #define ABI_FUNC 'qrfac'
 995 !End of the abilint section
 996 
 997  integer, intent(in)        :: m
 998  integer, intent(in)        :: n
 999  real (dp), intent(in out)  :: a(:,:)
1000  logical, intent(in)        :: pivot
1001  integer, intent(out)       :: ipvt(:)
1002  real (dp), intent(out)     :: rdiag(:)
1003  real (dp), intent(out)     :: acnorm(:)
1004 
1005  !  **********
1006 
1007  !  subroutine qrfac
1008 
1009  !  This subroutine uses Householder transformations with column pivoting
1010  !  (optional) to compute a qr factorization of the m by n matrix a.
1011  !  That is, qrfac determines an orthogonal matrix q, a permutation matrix p,
1012  !  and an upper trapezoidal matrix r with diagonal elements of nonincreasing
1013  !  magnitude, such that a*p = q*r.  The householder transformation for
1014  !  column k, k = 1,2,...,min(m,n), is of the form
1015 
1016  !                        t
1017  !        i - (1/u(k))*u*u
1018 
1019  !  where u has zeros in the first k-1 positions.  The form of this
1020  !  transformation and the method of pivoting first appeared in the
1021  !  corresponding linpack subroutine.
1022 
1023  !  the subroutine statement is
1024 
1025  !    subroutine qrfac(m, n, a, lda, pivot, ipvt, lipvt, rdiag, acnorm, wa)
1026 
1027  ! N.B. 3 of these arguments have been omitted in this version.
1028 
1029  !  where
1030 
1031  !    m is a positive integer input variable set to the number of rows of a.
1032 
1033  !    n is a positive integer input variable set to the number of columns of a.
1034 
1035  !    a is an m by n array.  On input a contains the matrix for
1036  !      which the qr factorization is to be computed.  On output
1037  !      the strict upper trapezoidal part of a contains the strict
1038  !      upper trapezoidal part of r, and the lower trapezoidal
1039  !      part of a contains a factored form of q (the non-trivial
1040  !      elements of the u vectors described above).
1041 
1042  !    lda is a positive integer input variable not less than m
1043  !      which specifies the leading dimension of the array a.
1044 
1045  !    pivot is a logical input variable.  If pivot is set true,
1046  !      then column pivoting is enforced.  If pivot is set false,
1047  !      then no column pivoting is done.
1048 
1049  !    ipvt is an integer output array of length lipvt.  ipvt
1050  !      defines the permutation matrix p such that a*p = q*r.
1051  !      Column j of p is column ipvt(j) of the identity matrix.
1052  !      If pivot is false, ipvt is not referenced.
1053 
1054  !    lipvt is a positive integer input variable.  If pivot is false,
1055  !      then lipvt may be as small as 1.  If pivot is true, then
1056  !      lipvt must be at least n.
1057 
1058  !    rdiag is an output array of length n which contains the
1059  !      diagonal elements of r.
1060 
1061  !    acnorm is an output array of length n which contains the norms of the
1062  !      corresponding columns of the input matrix a.
1063  !      If this information is not needed, then acnorm can coincide with rdiag.
1064 
1065  !    wa is a work array of length n.  If pivot is false, then wa
1066  !      can coincide with rdiag.
1067 
1068  !  subprograms called
1069 
1070  !    minpack-supplied ... dpmpar,enorm
1071 
1072  !    fortran-supplied ... MAX,SQRT,MIN
1073 
1074  !  argonne national laboratory. minpack project. march 1980.
1075  !  burton s. garbow, kenneth e. hillstrom, jorge j. more
1076 
1077  !  **********
1078  integer   :: i, j, jp1, k, kmax, minmn
1079  real (dp) :: ajnorm, epsmch, sum, temp, wa(n)
1080  real (dp), parameter :: one = 1.0_dp, p05 = 0.05_dp, zero = 0.0_dp
1081 
1082  !     epsmch is the machine precision.
1083 
1084  epsmch = EPSILON(zero)
1085 
1086  !     compute the initial column norms and initialize several arrays.
1087 
1088  DO  j = 1, n
1089    acnorm(j) = enorm(m,a(1:,j))
1090    rdiag(j) = acnorm(j)
1091    wa(j) = rdiag(j)
1092    IF (pivot) ipvt(j) = j
1093  END DO
1094 
1095  !     Reduce a to r with Householder transformations.
1096 
1097  minmn = MIN(m,n)
1098  DO  j = 1, minmn
1099    IF (.NOT.pivot) GO TO 40
1100 
1101  !        Bring the column of largest norm into the pivot position.
1102 
1103    kmax = j
1104    DO  k = j, n
1105      IF (rdiag(k) > rdiag(kmax)) kmax = k
1106    END DO
1107    IF (kmax == j) GO TO 40
1108    DO  i = 1, m
1109      temp = a(i,j)
1110      a(i,j) = a(i,kmax)
1111      a(i,kmax) = temp
1112    END DO
1113    rdiag(kmax) = rdiag(j)
1114    wa(kmax) = wa(j)
1115    k = ipvt(j)
1116    ipvt(j) = ipvt(kmax)
1117    ipvt(kmax) = k
1118 
1119  !     Compute the Householder transformation to reduce the
1120  !     j-th column of a to a multiple of the j-th unit vector.
1121 
1122    40 ajnorm = enorm(m-j+1, a(j:,j))
1123    IF (ajnorm == zero) CYCLE
1124    IF (a(j,j) < zero) ajnorm = -ajnorm
1125    a(j:m,j) = a(j:m,j)/ajnorm
1126    a(j,j) = a(j,j) + one
1127 
1128  !     Apply the transformation to the remaining columns and update the norms.
1129 
1130    jp1 = j + 1
1131    DO  k = jp1, n
1132      sum = DOT_PRODUCT( a(j:m,j), a(j:m,k) )
1133      temp = sum/a(j,j)
1134      a(j:m,k) = a(j:m,k) - temp*a(j:m,j)
1135      IF (.NOT.pivot .OR. rdiag(k) == zero) CYCLE
1136      temp = a(j,k)/rdiag(k)
1137      rdiag(k) = rdiag(k)*SQRT(MAX(zero, one-temp**2))
1138      IF (p05*(rdiag(k)/wa(k))**2 > epsmch) CYCLE
1139      rdiag(k) = enorm(m-j, a(jp1:,k))
1140      wa(k) = rdiag(k)
1141    END DO
1142    rdiag(j) = -ajnorm
1143  END DO
1144  RETURN
1145 
1146  !     last card of subroutine qrfac.
1147 
1148  end subroutine qrfac

m_levenberg_marquardt/qrsolv [ Functions ]

[ Top ] [ m_levenberg_marquardt ] [ Functions ]

NAME

FUNCTION

INPUTS

OUTPUT

NOTES

PARENTS

CHILDREN

SOURCE

1170  SUBROUTINE qrsolv(n, r, ipvt, diag, qtb, x, sdiag)
1171 
1172  ! N.B. Arguments LDR & WA have been removed.
1173 
1174 !This section has been created automatically by the script Abilint (TD).
1175 !Do not modify the following lines by hand.
1176 #undef ABI_FUNC
1177 #define ABI_FUNC 'qrsolv'
1178 !End of the abilint section
1179 
1180  integer, intent(in)        :: n
1181  real (dp), intent(in out)  :: r(:,:)
1182  integer, intent(in)        :: ipvt(:)
1183  real (dp), intent(in)      :: diag(:)
1184  real (dp), intent(in)      :: qtb(:)
1185  real (dp), intent(out)     :: x(:)
1186  real (dp), intent(out)     :: sdiag(:)
1187 
1188  !  **********
1189 
1190  !  subroutine qrsolv
1191 
1192  !  Given an m by n matrix a, an n by n diagonal matrix d, and an m-vector b,
1193  !  the problem is to determine an x which solves the system
1194 
1195  !        a*x = b ,     d*x = 0 ,
1196 
1197  !  in the least squares sense.
1198 
1199  !  This subroutine completes the solution of the problem if it is provided
1200  !  with the necessary information from the qr factorization, with column
1201  !  pivoting, of a.  That is, if a*p = q*r, where p is a permutation matrix,
1202  !  q has orthogonal columns, and r is an upper triangular matrix with diagonal
1203  !  elements of nonincreasing magnitude, then qrsolv expects the full upper
1204  !  triangle of r, the permutation matrix p, and the first n components of
1205  !  (q transpose)*b.  The system a*x = b, d*x = 0, is then equivalent to
1206 
1207  !               t       t
1208  !        r*z = q *b ,  p *d*p*z = 0 ,
1209 
1210  !  where x = p*z. if this system does not have full rank,
1211  !  then a least squares solution is obtained.  On output qrsolv
1212  !  also provides an upper triangular matrix s such that
1213 
1214  !         t   t               t
1215  !        p *(a *a + d*d)*p = s *s .
1216 
1217  !  s is computed within qrsolv and may be of separate interest.
1218 
1219  !  the subroutine statement is
1220 
1221  !    subroutine qrsolv(n, r, ldr, ipvt, diag, qtb, x, sdiag, wa)
1222 
1223  ! N.B. Arguments LDR and WA have been removed in this version.
1224 
1225  !  where
1226 
1227  !    n is a positive integer input variable set to the order of r.
1228 
1229  !    r is an n by n array.  On input the full upper triangle must contain
1230  !      the full upper triangle of the matrix r.
1231  !      On output the full upper triangle is unaltered, and the strict lower
1232  !      triangle contains the strict upper triangle (transposed) of the
1233  !      upper triangular matrix s.
1234 
1235  !    ldr is a positive integer input variable not less than n
1236  !      which specifies the leading dimension of the array r.
1237 
1238  !    ipvt is an integer input array of length n which defines the
1239  !      permutation matrix p such that a*p = q*r.  Column j of p
1240  !      is column ipvt(j) of the identity matrix.
1241 
1242  !    diag is an input array of length n which must contain the
1243  !      diagonal elements of the matrix d.
1244 
1245  !    qtb is an input array of length n which must contain the first
1246  !      n elements of the vector (q transpose)*b.
1247 
1248  !    x is an output array of length n which contains the least
1249  !      squares solution of the system a*x = b, d*x = 0.
1250 
1251  !    sdiag is an output array of length n which contains the
1252  !      diagonal elements of the upper triangular matrix s.
1253 
1254  !    wa is a work array of length n.
1255 
1256  !  subprograms called
1257 
1258  !    fortran-supplied ... ABS,SQRT
1259 
1260  !  argonne national laboratory. minpack project. march 1980.
1261  !  burton s. garbow, kenneth e. hillstrom, jorge j. more
1262 
1263  !  **********
1264  integer   :: i, j, k, kp1, l, nsing
1265  real (dp) :: cos, cotan, qtbpj, sin, sum, tan, temp, wa(n)
1266  real (dp), parameter :: p5 = 0.5_dp, p25 = 0.25_dp, zero = 0.0_dp
1267 
1268  !     Copy r and (q transpose)*b to preserve input and initialize s.
1269  !     In particular, save the diagonal elements of r in x.
1270 
1271  DO  j = 1, n
1272    r(j:n,j) = r(j,j:n)
1273    x(j) = r(j,j)
1274    wa(j) = qtb(j)
1275  END DO
1276 
1277  !     Eliminate the diagonal matrix d using a givens rotation.
1278 
1279  DO  j = 1, n
1280 
1281  !        Prepare the row of d to be eliminated, locating the
1282  !        diagonal element using p from the qr factorization.
1283 
1284    l = ipvt(j)
1285    IF (diag(l) == zero) CYCLE
1286    sdiag(j:n) = zero
1287    sdiag(j) = diag(l)
1288 
1289  !     The transformations to eliminate the row of d modify only a single
1290  !     element of (q transpose)*b beyond the first n, which is initially zero.
1291 
1292    qtbpj = zero
1293    DO  k = j, n
1294 
1295  !        Determine a givens rotation which eliminates the
1296  !        appropriate element in the current row of d.
1297 
1298      IF (sdiag(k) == zero) CYCLE
1299      IF (ABS(r(k,k)) < ABS(sdiag(k))) THEN
1300        cotan = r(k,k)/sdiag(k)
1301        SIN = p5/SQRT(p25 + p25*cotan**2)
1302        COS = SIN*cotan
1303      ELSE
1304        TAN = sdiag(k)/r(k,k)
1305        COS = p5/SQRT(p25 + p25*TAN**2)
1306        SIN = COS*TAN
1307      END IF
1308 
1309  !        Compute the modified diagonal element of r and
1310  !        the modified element of ((q transpose)*b,0).
1311 
1312      r(k,k) = COS*r(k,k) + SIN*sdiag(k)
1313      temp = COS*wa(k) + SIN*qtbpj
1314      qtbpj = -SIN*wa(k) + COS*qtbpj
1315      wa(k) = temp
1316 
1317  !        Accumulate the tranformation in the row of s.
1318 
1319      kp1 = k + 1
1320      DO  i = kp1, n
1321        temp = COS*r(i,k) + SIN*sdiag(i)
1322        sdiag(i) = -SIN*r(i,k) + COS*sdiag(i)
1323        r(i,k) = temp
1324      END DO
1325    END DO
1326 
1327  !     Store the diagonal element of s and restore
1328  !     the corresponding diagonal element of r.
1329 
1330    sdiag(j) = r(j,j)
1331    r(j,j) = x(j)
1332  END DO
1333 
1334  !     Solve the triangular system for z.  If the system is singular,
1335  !     then obtain a least squares solution.
1336 
1337  nsing = n
1338  DO  j = 1, n
1339    IF (sdiag(j) == zero .AND. nsing == n) nsing = j - 1
1340    IF (nsing < n) wa(j) = zero
1341  END DO
1342 
1343  DO  k = 1, nsing
1344    j = nsing - k + 1
1345    sum = DOT_PRODUCT( r(j+1:nsing,j), wa(j+1:nsing) )
1346    wa(j) = (wa(j) - sum)/sdiag(j)
1347  END DO
1348 
1349  !     Permute the components of z back to components of x.
1350 
1351  DO  j = 1, n
1352    l = ipvt(j)
1353    x(l) = wa(j)
1354  END DO
1355  RETURN
1356 
1357  !     last card of subroutine qrsolv.
1358 
1359  END SUBROUTINE qrsolv