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