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) 8c ********** 9c 10c subroutine lmpar 11c 12c given an m by n matrix a, an n by n nonsingular diagonal 13c matrix d, an m-vector b, and a positive number delta, 14c the problem is to determine a value for the parameter 15c par such that if x solves the system 16c 17c a*x = b , sqrt(par)*d*x = 0 , 18c 19c in the least squares sense, and dxnorm is the euclidean 20c norm of d*x, then either par is zero and 21c 22c (dxnorm-delta) .le. 0.1*delta , 23c 24c or par is positive and 25c 26c abs(dxnorm-delta) .le. 0.1*delta . 27c 28c this subroutine completes the solution of the problem 29c if it is provided with the necessary information from the 30c qr factorization, with column pivoting, of a. that is, if 31c a*p = q*r, where p is a permutation matrix, q has orthogonal 32c columns, and r is an upper triangular matrix with diagonal 33c elements of nonincreasing magnitude, then lmpar expects 34c the full upper triangle of r, the permutation matrix p, 35c and the first n components of (q transpose)*b. on output 36c lmpar also provides an upper triangular matrix s such that 37c 38c t t t 39c p *(a *a + par*d*d)*p = s *s . 40c 41c s is employed within lmpar and may be of separate interest. 42c 43c only a few iterations are generally needed for convergence 44c of the algorithm. if, however, the limit of 10 iterations 45c is reached, then the output par will contain the best 46c value obtained so far. 47c 48c the subroutine statement is 49c 50c subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag, 51c wa1,wa2) 52c 53c where 54c 55c n is a positive integer input variable set to the order of r. 56c 57c r is an n by n array. on input the full upper triangle 58c must contain the full upper triangle of the matrix r. 59c on output the full upper triangle is unaltered, and the 60c strict lower triangle contains the strict upper triangle 61c (transposed) of the upper triangular matrix s. 62c 63c ldr is a positive integer input variable not less than n 64c which specifies the leading dimension of the array r. 65c 66c ipvt is an integer input array of length n which defines the 67c permutation matrix p such that a*p = q*r. column j of p 68c is column ipvt(j) of the identity matrix. 69c 70c diag is an input array of length n which must contain the 71c diagonal elements of the matrix d. 72c 73c qtb is an input array of length n which must contain the first 74c n elements of the vector (q transpose)*b. 75c 76c delta is a positive input variable which specifies an upper 77c bound on the euclidean norm of d*x. 78c 79c par is a nonnegative variable. on input par contains an 80c initial estimate of the levenberg-marquardt parameter. 81c on output par contains the final estimate. 82c 83c x is an output array of length n which contains the least 84c squares solution of the system a*x = b, sqrt(par)*d*x = 0, 85c for the output par. 86c 87c sdiag is an output array of length n which contains the 88c diagonal elements of the upper triangular matrix s. 89c 90c wa1 and wa2 are work arrays of length n. 91c 92c subprograms called 93c 94c minpack-supplied ... dpmpar,enorm,qrsolv 95c 96c fortran-supplied ... dabs,dmax1,dmin1,dsqrt 97c 98c argonne national laboratory. minpack project. march 1980. 99c burton s. garbow, kenneth e. hillstrom, jorge j. more 100c 101c ********** 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/ 107c 108c dwarf is the smallest positive magnitude. 109c 110 dwarf = dpmpar(2) 111c 112c compute and store in x the gauss-newton direction. if the 113c jacobian is rank-deficient, obtain a least squares solution. 114c 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 138c 139c initialize the iteration counter. 140c evaluate the function at the origin, and test 141c for acceptance of the gauss-newton direction. 142c 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 150c 151c if the jacobian is not rank deficient, the newton 152c step provides a lower bound, parl, for the zero of 153c the function. otherwise set this bound to zero. 154c 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 174c 175c calculate an upper bound, paru, for the zero of the function. 176c 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) 188c 189c if the input par lies outside of the interval (parl,paru), 190c set par to the closer endpoint. 191c 192 par = dmax1(par,parl) 193 par = dmin1(par,paru) 194 if (par .eq. zero) par = gnorm/dxnorm 195c 196c beginning of an iteration. 197c 198 150 continue 199 iter = iter + 1 200c 201c evaluate the function at the current value of par. 202c 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 215c 216c if the function is small enough, accept the current value 217c of par. also test for the exceptional cases where parl 218c is zero or the number of iterations has reached 10. 219c 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 223c 224c compute the newton correction. 225c 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 242c 243c depending on the sign of the function, update parl or paru. 244c 245 if (fp .gt. zero) parl = dmax1(parl,par) 246 if (fp .lt. zero) paru = dmin1(paru,par) 247c 248c compute an improved estimate for par. 249c 250 par = dmax1(parl,par+parc) 251c 252c end of an iteration. 253c 254 go to 150 255 220 continue 256c 257c termination. 258c 259 if (iter .eq. 0) par = zero 260 return 261c 262c last card of subroutine lmpar. 263c 264 end 265