share/minpack/fortran/dogleg.f

   1       subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2)
   2       integer n,lr
   3       double precision delta
   4       double precision r(lr),diag(n),qtb(n),x(n),wa1(n),wa2(n)
   5 c     **********
   6 c
   7 c     subroutine dogleg
   8 c
   9 c     given an m by n matrix a, an n by n nonsingular diagonal
  10 c     matrix d, an m-vector b, and a positive number delta, the
  11 c     problem is to determine the convex combination x of the
  12 c     gauss-newton and scaled gradient directions that minimizes
  13 c     (a*x - b) in the least squares sense, subject to the
  14 c     restriction that the euclidean norm of d*x be at most delta.
  15 c
  16 c     this subroutine completes the solution of the problem
  17 c     if it is provided with the necessary information from the
  18 c     qr factorization of a. that is, if a = q*r, where q has
  19 c     orthogonal columns and r is an upper triangular matrix,
  20 c     then dogleg expects the full upper triangle of r and
  21 c     the first n components of (q transpose)*b.
  22 c
  23 c     the subroutine statement is
  24 c
  25 c       subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2)
  26 c
  27 c     where
  28 c
  29 c       n is a positive integer input variable set to the order of r.
  30 c
  31 c       r is an input array of length lr which must contain the upper
  32 c         triangular matrix r stored by rows.
  33 c
  34 c       lr is a positive integer input variable not less than
  35 c         (n*(n+1))/2.
  36 c
  37 c       diag is an input array of length n which must contain the
  38 c         diagonal elements of the matrix d.
  39 c
  40 c       qtb is an input array of length n which must contain the first
  41 c         n elements of the vector (q transpose)*b.
  42 c
  43 c       delta is a positive input variable which specifies an upper
  44 c         bound on the euclidean norm of d*x.
  45 c
  46 c       x is an output array of length n which contains the desired
  47 c         convex combination of the gauss-newton direction and the
  48 c         scaled gradient direction.
  49 c
  50 c       wa1 and wa2 are work arrays of length n.
  51 c
  52 c     subprograms called
  53 c
  54 c       minpack-supplied ... dpmpar,enorm
  55 c
  56 c       fortran-supplied ... dabs,dmax1,dmin1,dsqrt
  57 c
  58 c     argonne national laboratory. minpack project. march 1980.
  59 c     burton s. garbow, kenneth e. hillstrom, jorge j. more
  60 c
  61 c     **********
  62       integer i,j,jj,jp1,k,l
  63       double precision alpha,bnorm,epsmch,gnorm,one,qnorm,sgnorm,sum,
  64      *                 temp,zero
  65       double precision dpmpar,enorm
  66       data one,zero /1.0d0,0.0d0/
  67 c
  68 c     epsmch is the machine precision.
  69 c
  70       epsmch = dpmpar(1)
  71 c
  72 c     first, calculate the gauss-newton direction.
  73 c
  74       jj = (n*(n + 1))/2 + 1
  75       do 50 k = 1, n
  76          j = n - k + 1
  77          jp1 = j + 1
  78          jj = jj - k
  79          l = jj + 1
  80          sum = zero
  81          if (n .lt. jp1) go to 20
  82          do 10 i = jp1, n
  83             sum = sum + r(l)*x(i)
  84             l = l + 1
  85    10       continue
  86    20    continue
  87          temp = r(jj)
  88          if (temp .ne. zero) go to 40
  89          l = j
  90          do 30 i = 1, j
  91             temp = dmax1(temp,dabs(r(l)))
  92             l = l + n - i
  93    30       continue
  94          temp = epsmch*temp
  95          if (temp .eq. zero) temp = epsmch
  96    40    continue
  97          x(j) = (qtb(j) - sum)/temp
  98    50    continue
  99 c
 100 c     test whether the gauss-newton direction is acceptable.
 101 c
 102       do 60 j = 1, n
 103          wa1(j) = zero
 104          wa2(j) = diag(j)*x(j)
 105    60    continue
 106       qnorm = enorm(n,wa2)
 107       if (qnorm .le. delta) go to 140
 108 c
 109 c     the gauss-newton direction is not acceptable.
 110 c     next, calculate the scaled gradient direction.
 111 c
 112       l = 1
 113       do 80 j = 1, n
 114          temp = qtb(j)
 115          do 70 i = j, n
 116             wa1(i) = wa1(i) + r(l)*temp
 117             l = l + 1
 118    70       continue
 119          wa1(j) = wa1(j)/diag(j)
 120    80    continue
 121 c
 122 c     calculate the norm of the scaled gradient and test for
 123 c     the special case in which the scaled gradient is zero.
 124 c
 125       gnorm = enorm(n,wa1)
 126       sgnorm = zero
 127       alpha = delta/qnorm
 128       if (gnorm .eq. zero) go to 120
 129 c
 130 c     calculate the point along the scaled gradient
 131 c     at which the quadratic is minimized.
 132 c
 133       do 90 j = 1, n
 134          wa1(j) = (wa1(j)/gnorm)/diag(j)
 135    90    continue
 136       l = 1
 137       do 110 j = 1, n
 138          sum = zero
 139          do 100 i = j, n
 140             sum = sum + r(l)*wa1(i)
 141             l = l + 1
 142   100       continue
 143          wa2(j) = sum
 144   110    continue
 145       temp = enorm(n,wa2)
 146       sgnorm = (gnorm/temp)/temp
 147 c
 148 c     test whether the scaled gradient direction is acceptable.
 149 c
 150       alpha = zero
 151       if (sgnorm .ge. delta) go to 120
 152 c
 153 c     the scaled gradient direction is not acceptable.
 154 c     finally, calculate the point along the dogleg
 155 c     at which the quadratic is minimized.
 156 c
 157       bnorm = enorm(n,qtb)
 158       temp = (bnorm/gnorm)*(bnorm/qnorm)*(sgnorm/delta)
 159       temp = temp - (delta/qnorm)*(sgnorm/delta)**2
 160      *       + dsqrt((temp-(delta/qnorm))**2
 161      *               +(one-(delta/qnorm)**2)*(one-(sgnorm/delta)**2))
 162       alpha = ((delta/qnorm)*(one - (sgnorm/delta)**2))/temp
 163   120 continue
 164 c
 165 c     form appropriate convex combination of the gauss-newton
 166 c     direction and the scaled gradient direction.
 167 c
 168       temp = (one - alpha)*dmin1(sgnorm,delta)
 169       do 130 j = 1, n
 170          x(j) = temp*wa1(j) + alpha*x(j)
 171   130    continue
 172   140 continue
 173       return
 174 c
 175 c     last card of subroutine dogleg.
 176 c
 177       end