Mercurial > hg > octave-nkf

comparison scripts/optimization/fsolve.m @ 8990:349d75161672

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
more cosmetic adjustments to fsolve
author Jaroslav Hajek <highegg@gmail.com>
date Tue, 17 Mar 2009 13:38:16 +0100
parents 315828058e0d
children c235a59d30a4
comparison
equal deleted inserted replaced
8989:46a12e3f882c 8990:349d75161672
162 162
163 tolx = optimget (options, "TolX", sqrt (macheps)); 163 tolx = optimget (options, "TolX", sqrt (macheps));
164 tolf = optimget (options, "TolFun", sqrt (macheps)); 164 tolf = optimget (options, "TolFun", sqrt (macheps));
165 165
166 factor = 100; 166 factor = 100;
167 ## FIXME: TypicalX corresponds to user scaling (???)
168 autodg = true;
169 167
170 niter = 1; 168 niter = 1;
171 nfev = 1; 169 nfev = 1;
172 170
173 x = x0(:); 171 x = x0(:);
174 info = 0; 172 info = 0;
175 173
176 ## Initial evaluation. 174 ## Initial evaluation.
175 ## Handle arbitrary shapes of x and f and remember them.
177 fvec = fcn (reshape (x, xsiz)); 176 fvec = fcn (reshape (x, xsiz));
178 fsiz = size (fvec); 177 fsiz = size (fvec);
179 fvec = fvec(:); 178 fvec = fvec(:);
180 fn = norm (fvec); 179 fn = norm (fvec);
181 m = length (fvec); 180 m = length (fvec);
198 197
199 ## Outer loop. 198 ## Outer loop.
200 while (niter < maxiter && nfev < maxfev && ! info) 199 while (niter < maxiter && nfev < maxfev && ! info)
201 200
202 ## Calculate function value and Jacobian (possibly via FD). 201 ## Calculate function value and Jacobian (possibly via FD).
203 ## Handle arbitrary shapes of x and f and remember them.
204 if (has_jac) 202 if (has_jac)
205 [fvec, fjac] = fcn (reshape (x, xsiz)); 203 [fvec, fjac] = fcn (reshape (x, xsiz));
206 ## If the jacobian is sparse, disable Broyden updating. 204 ## If the jacobian is sparse, disable Broyden updating.
207 if (issparse (fjac)) 205 if (issparse (fjac))
208 updating = false; 206 updating = false;
229 endif 227 endif
230 228
231 ## Get column norms, use them as scaling factors. 229 ## Get column norms, use them as scaling factors.
232 jcn = norm (fjac, 'columns').'; 230 jcn = norm (fjac, 'columns').';
233 if (niter == 1) 231 if (niter == 1)
234 if (autodg) 232 dg = jcn;
235 dg = jcn; 233 dg(dg == 0) = 1;
236 dg(dg == 0) = 1;
237 endif
238 xn = norm (dg .* x); 234 xn = norm (dg .* x);
239 ## FIXME: something better? 235 ## FIXME: something better?
240 delta = max (factor * xn, 1); 236 delta = max (factor * xn, 1);
241 endif 237 endif
242 238
243 ## Rescale if necessary. 239 ## Rescale adaptively.
244 if (autodg) 240 ## FIXME: the original minpack used the following rescaling strategy:
245 ## FIXME: the original minpack used the following rescaling strategy: 241 ## dg = max (dg, jcn);
246 ## dg = max (dg, jcn); 242 ## but it seems not good if we start with a bad guess yielding jacobian
247 ## but it seems not good if we start with a bad guess yielding jacobian 243 ## columns with large norms that later decrease, because the corresponding
248 ## columns with large norms that later decrease, because the corresponding 244 ## variable will still be overscaled. So instead, we only give the old
249 ## variable will still be overscaled. So instead, we only give the old 245 ## scaling a small momentum, but do not honor it.
250 ## scaling a small momentum, but do not honor it. 246
251 247 dg = max (0.1*dg, jcn);
252 dg = max (0.1*dg, jcn); 248
253 249 ## It also seems that in the case of fast (and inhomogeneously) changing
254 ## It also seems that in the case of fast (and inhomogeneously) changing 250 ## jacobian, the Broyden updates are of little use, so maybe we could
255 ## jacobian, the Broyden updates are of little use, so maybe we could 251 ## skip them if a big disproportional change is expected. The question is,
256 ## skip them if a big disproportional change is expected. The question is, 252 ## of course, how to define the above terms :)
257 ## of course, how to define the above terms :)
258 endif
259 253
260 lastratio = 0; 254 lastratio = 0;
261 nfail = 0; 255 nfail = 0;
262 nsuc = 0; 256 nsuc = 0;
263 ## Inner loop. 257 ## Inner loop.