share/minpack/fortran/lmpar.f

   1       subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag,wa1,
   2      *                 wa2)
   3       integer n,ldr
   4       integer ipvt(n)
   5       double precision delta,par
   6       double precision r(ldr,n),diag(n),qtb(n),x(n),sdiag(n),wa1(n),
   7      *                 wa2(n)
   8 c     **********
   9 c
  10 c     subroutine lmpar
  11 c
  12 c     given an m by n matrix a, an n by n nonsingular diagonal
  13 c     matrix d, an m-vector b, and a positive number delta,
  14 c     the problem is to determine a value for the parameter
  15 c     par such that if x solves the system
  16 c
  17 c           a*x = b ,     sqrt(par)*d*x = 0 ,
  18 c
  19 c     in the least squares sense, and dxnorm is the euclidean
  20 c     norm of d*x, then either par is zero and
  21 c
  22 c           (dxnorm-delta) .le. 0.1*delta ,
  23 c
  24 c     or par is positive and
  25 c
  26 c           abs(dxnorm-delta) .le. 0.1*delta .
  27 c
  28 c     this subroutine completes the solution of the problem
  29 c     if it is provided with the necessary information from the
  30 c     qr factorization, with column pivoting, of a. that is, if
  31 c     a*p = q*r, where p is a permutation matrix, q has orthogonal
  32 c     columns, and r is an upper triangular matrix with diagonal
  33 c     elements of nonincreasing magnitude, then lmpar expects
  34 c     the full upper triangle of r, the permutation matrix p,
  35 c     and the first n components of (q transpose)*b. on output
  36 c     lmpar also provides an upper triangular matrix s such that
  37 c
  38 c            t   t                   t
  39 c           p *(a *a + par*d*d)*p = s *s .
  40 c
  41 c     s is employed within lmpar and may be of separate interest.
  42 c
  43 c     only a few iterations are generally needed for convergence
  44 c     of the algorithm. if, however, the limit of 10 iterations
  45 c     is reached, then the output par will contain the best
  46 c     value obtained so far.
  47 c
  48 c     the subroutine statement is
  49 c
  50 c       subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag,
  51 c                        wa1,wa2)
  52 c
  53 c     where
  54 c
  55 c       n is a positive integer input variable set to the order of r.
  56 c
  57 c       r is an n by n array. on input the full upper triangle
  58 c         must contain the full upper triangle of the matrix r.
  59 c         on output the full upper triangle is unaltered, and the
  60 c         strict lower triangle contains the strict upper triangle
  61 c         (transposed) of the upper triangular matrix s.
  62 c
  63 c       ldr is a positive integer input variable not less than n
  64 c         which specifies the leading dimension of the array r.
  65 c
  66 c       ipvt is an integer input array of length n which defines the
  67 c         permutation matrix p such that a*p = q*r. column j of p
  68 c         is column ipvt(j) of the identity matrix.
  69 c
  70 c       diag is an input array of length n which must contain the
  71 c         diagonal elements of the matrix d.
  72 c
  73 c       qtb is an input array of length n which must contain the first
  74 c         n elements of the vector (q transpose)*b.
  75 c
  76 c       delta is a positive input variable which specifies an upper
  77 c         bound on the euclidean norm of d*x.
  78 c
  79 c       par is a nonnegative variable. on input par contains an
  80 c         initial estimate of the levenberg-marquardt parameter.
  81 c         on output par contains the final estimate.
  82 c
  83 c       x is an output array of length n which contains the least
  84 c         squares solution of the system a*x = b, sqrt(par)*d*x = 0,
  85 c         for the output par.
  86 c
  87 c       sdiag is an output array of length n which contains the
  88 c         diagonal elements of the upper triangular matrix s.
  89 c
  90 c       wa1 and wa2 are work arrays of length n.
  91 c
  92 c     subprograms called
  93 c
  94 c       minpack-supplied ... dpmpar,enorm,qrsolv
  95 c
  96 c       fortran-supplied ... dabs,dmax1,dmin1,dsqrt
  97 c
  98 c     argonne national laboratory. minpack project. march 1980.
  99 c     burton s. garbow, kenneth e. hillstrom, jorge j. more
 100 c
 101 c     **********
 102       integer i,iter,j,jm1,jp1,k,l,nsing
 103       double precision dxnorm,dwarf,fp,gnorm,parc,parl,paru,p1,p001,
 104      *                 sum,temp,zero
 105       double precision dpmpar,enorm
 106       data p1,p001,zero /1.0d-1,1.0d-3,0.0d0/
 107 c
 108 c     dwarf is the smallest positive magnitude.
 109 c
 110       dwarf = dpmpar(2)
 111 c
 112 c     compute and store in x the gauss-newton direction. if the
 113 c     jacobian is rank-deficient, obtain a least squares solution.
 114 c
 115       nsing = n
 116       do 10 j = 1, n
 117          wa1(j) = qtb(j)
 118          if (r(j,j) .eq. zero .and. nsing .eq. n) nsing = j - 1
 119          if (nsing .lt. n) wa1(j) = zero
 120    10    continue
 121       if (nsing .lt. 1) go to 50
 122       do 40 k = 1, nsing
 123          j = nsing - k + 1
 124          wa1(j) = wa1(j)/r(j,j)
 125          temp = wa1(j)
 126          jm1 = j - 1
 127          if (jm1 .lt. 1) go to 30
 128          do 20 i = 1, jm1
 129             wa1(i) = wa1(i) - r(i,j)*temp
 130    20       continue
 131    30    continue
 132    40    continue
 133    50 continue
 134       do 60 j = 1, n
 135          l = ipvt(j)
 136          x(l) = wa1(j)
 137    60    continue
 138 c
 139 c     initialize the iteration counter.
 140 c     evaluate the function at the origin, and test
 141 c     for acceptance of the gauss-newton direction.
 142 c
 143       iter = 0
 144       do 70 j = 1, n
 145          wa2(j) = diag(j)*x(j)
 146    70    continue
 147       dxnorm = enorm(n,wa2)
 148       fp = dxnorm - delta
 149       if (fp .le. p1*delta) go to 220
 150 c
 151 c     if the jacobian is not rank deficient, the newton
 152 c     step provides a lower bound, parl, for the zero of
 153 c     the function. otherwise set this bound to zero.
 154 c
 155       parl = zero
 156       if (nsing .lt. n) go to 120
 157       do 80 j = 1, n
 158          l = ipvt(j)
 159          wa1(j) = diag(l)*(wa2(l)/dxnorm)
 160    80    continue
 161       do 110 j = 1, n
 162          sum = zero
 163          jm1 = j - 1
 164          if (jm1 .lt. 1) go to 100
 165          do 90 i = 1, jm1
 166             sum = sum + r(i,j)*wa1(i)
 167    90       continue
 168   100    continue
 169          wa1(j) = (wa1(j) - sum)/r(j,j)
 170   110    continue
 171       temp = enorm(n,wa1)
 172       parl = ((fp/delta)/temp)/temp
 173   120 continue
 174 c
 175 c     calculate an upper bound, paru, for the zero of the function.
 176 c
 177       do 140 j = 1, n
 178          sum = zero
 179          do 130 i = 1, j
 180             sum = sum + r(i,j)*qtb(i)
 181   130       continue
 182          l = ipvt(j)
 183          wa1(j) = sum/diag(l)
 184   140    continue
 185       gnorm = enorm(n,wa1)
 186       paru = gnorm/delta
 187       if (paru .eq. zero) paru = dwarf/dmin1(delta,p1)
 188 c
 189 c     if the input par lies outside of the interval (parl,paru),
 190 c     set par to the closer endpoint.
 191 c
 192       par = dmax1(par,parl)
 193       par = dmin1(par,paru)
 194       if (par .eq. zero) par = gnorm/dxnorm
 195 c
 196 c     beginning of an iteration.
 197 c
 198   150 continue
 199          iter = iter + 1
 200 c
 201 c        evaluate the function at the current value of par.
 202 c
 203          if (par .eq. zero) par = dmax1(dwarf,p001*paru)
 204          temp = dsqrt(par)
 205          do 160 j = 1, n
 206             wa1(j) = temp*diag(j)
 207   160       continue
 208          call qrsolv(n,r,ldr,ipvt,wa1,qtb,x,sdiag,wa2)
 209          do 170 j = 1, n
 210             wa2(j) = diag(j)*x(j)
 211   170       continue
 212          dxnorm = enorm(n,wa2)
 213          temp = fp
 214          fp = dxnorm - delta
 215 c
 216 c        if the function is small enough, accept the current value
 217 c        of par. also test for the exceptional cases where parl
 218 c        is zero or the number of iterations has reached 10.
 219 c
 220          if (dabs(fp) .le. p1*delta
 221      *       .or. parl .eq. zero .and. fp .le. temp
 222      *            .and. temp .lt. zero .or. iter .eq. 10) go to 220
 223 c
 224 c        compute the newton correction.
 225 c
 226          do 180 j = 1, n
 227             l = ipvt(j)
 228             wa1(j) = diag(l)*(wa2(l)/dxnorm)
 229   180       continue
 230          do 210 j = 1, n
 231             wa1(j) = wa1(j)/sdiag(j)
 232             temp = wa1(j)
 233             jp1 = j + 1
 234             if (n .lt. jp1) go to 200
 235             do 190 i = jp1, n
 236                wa1(i) = wa1(i) - r(i,j)*temp
 237   190          continue
 238   200       continue
 239   210       continue
 240          temp = enorm(n,wa1)
 241          parc = ((fp/delta)/temp)/temp
 242 c
 243 c        depending on the sign of the function, update parl or paru.
 244 c
 245          if (fp .gt. zero) parl = dmax1(parl,par)
 246          if (fp .lt. zero) paru = dmin1(paru,par)
 247 c
 248 c        compute an improved estimate for par.
 249 c
 250          par = dmax1(parl,par+parc)
 251 c
 252 c        end of an iteration.
 253 c
 254          go to 150
 255   220 continue
 256 c
 257 c     termination.
 258 c
 259       if (iter .eq. 0) par = zero
 260       return
 261 c
 262 c     last card of subroutine lmpar.
 263 c
 264       end