Line data Source code
1 0 : subroutine pda_lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag,wa1,
2 0 : * 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 pda_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 pda_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 pda_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 pda_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 pda_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 ... pda_dpmpar,pda_enorm,pda_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 pda_dpmpar,pda_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 0 : dwarf = pda_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 0 : nsing = n
116 0 : do 10 j = 1, n
117 0 : wa1(j) = qtb(j)
118 0 : if (r(j,j) .eq. zero .and. nsing .eq. n) nsing = j - 1
119 0 : if (nsing .lt. n) wa1(j) = zero
120 0 : 10 continue
121 0 : if (nsing .lt. 1) go to 50
122 0 : do 40 k = 1, nsing
123 0 : j = nsing - k + 1
124 0 : wa1(j) = wa1(j)/r(j,j)
125 0 : temp = wa1(j)
126 0 : jm1 = j - 1
127 0 : if (jm1 .lt. 1) go to 30
128 0 : do 20 i = 1, jm1
129 0 : wa1(i) = wa1(i) - r(i,j)*temp
130 0 : 20 continue
131 : 30 continue
132 0 : 40 continue
133 : 50 continue
134 0 : do 60 j = 1, n
135 0 : l = ipvt(j)
136 0 : x(l) = wa1(j)
137 0 : 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 0 : iter = 0
144 0 : do 70 j = 1, n
145 0 : wa2(j) = diag(j)*x(j)
146 0 : 70 continue
147 0 : dxnorm = pda_enorm(n,wa2)
148 0 : fp = dxnorm - delta
149 0 : 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 0 : parl = zero
156 0 : if (nsing .lt. n) go to 120
157 0 : do 80 j = 1, n
158 0 : l = ipvt(j)
159 0 : wa1(j) = diag(l)*(wa2(l)/dxnorm)
160 0 : 80 continue
161 0 : do 110 j = 1, n
162 0 : sum = zero
163 0 : jm1 = j - 1
164 0 : if (jm1 .lt. 1) go to 100
165 0 : do 90 i = 1, jm1
166 0 : sum = sum + r(i,j)*wa1(i)
167 0 : 90 continue
168 : 100 continue
169 0 : wa1(j) = (wa1(j) - sum)/r(j,j)
170 0 : 110 continue
171 0 : temp = pda_enorm(n,wa1)
172 0 : 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 0 : do 140 j = 1, n
178 0 : sum = zero
179 0 : do 130 i = 1, j
180 0 : sum = sum + r(i,j)*qtb(i)
181 0 : 130 continue
182 0 : l = ipvt(j)
183 0 : wa1(j) = sum/diag(l)
184 0 : 140 continue
185 0 : gnorm = pda_enorm(n,wa1)
186 0 : paru = gnorm/delta
187 0 : 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 0 : par = dmax1(par,parl)
193 0 : par = dmin1(par,paru)
194 0 : if (par .eq. zero) par = gnorm/dxnorm
195 : c
196 : c beginning of an iteration.
197 : c
198 : 150 continue
199 0 : iter = iter + 1
200 : c
201 : c evaluate the function at the current value of par.
202 : c
203 0 : if (par .eq. zero) par = dmax1(dwarf,p001*paru)
204 0 : temp = dsqrt(par)
205 0 : do 160 j = 1, n
206 0 : wa1(j) = temp*diag(j)
207 0 : 160 continue
208 0 : call pda_qrsolv(n,r,ldr,ipvt,wa1,qtb,x,sdiag,wa2)
209 0 : do 170 j = 1, n
210 0 : wa2(j) = diag(j)*x(j)
211 0 : 170 continue
212 0 : dxnorm = pda_enorm(n,wa2)
213 0 : temp = fp
214 0 : 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 0 : * .and. temp .lt. zero .or. iter .eq. 10) go to 220
223 : c
224 : c compute the newton correction.
225 : c
226 0 : do 180 j = 1, n
227 0 : l = ipvt(j)
228 0 : wa1(j) = diag(l)*(wa2(l)/dxnorm)
229 0 : 180 continue
230 0 : do 210 j = 1, n
231 0 : wa1(j) = wa1(j)/sdiag(j)
232 0 : temp = wa1(j)
233 0 : jp1 = j + 1
234 0 : if (n .lt. jp1) go to 200
235 0 : do 190 i = jp1, n
236 0 : wa1(i) = wa1(i) - r(i,j)*temp
237 0 : 190 continue
238 : 200 continue
239 0 : 210 continue
240 0 : temp = pda_enorm(n,wa1)
241 0 : parc = ((fp/delta)/temp)/temp
242 : c
243 : c depending on the sign of the function, update parl or paru.
244 : c
245 0 : if (fp .gt. zero) parl = dmax1(parl,par)
246 0 : if (fp .lt. zero) paru = dmin1(paru,par)
247 : c
248 : c compute an improved estimate for par.
249 : c
250 0 : par = dmax1(parl,par+parc)
251 : c
252 : c end of an iteration.
253 : c
254 0 : go to 150
255 : 220 continue
256 : c
257 : c termination.
258 : c
259 0 : if (iter .eq. 0) par = zero
260 0 : return
261 : c
262 : c last card of subroutine pda_lmpar.
263 : c
264 : end
|
