← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to scipy/optimize/minpack/hybrd.f

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
scipy/optimize/minpack/hybrd.f
 
1
      subroutine hybrd(fcn,n,x,fvec,xtol,maxfev,ml,mu,epsfcn,diag,
 
2
     *                 mode,factor,nprint,info,nfev,fjac,ldfjac,r,lr,
 
3
     *                 qtf,wa1,wa2,wa3,wa4)
 
4
      integer n,maxfev,ml,mu,mode,nprint,info,nfev,ldfjac,lr
 
5
      double precision xtol,epsfcn,factor
 
6
      double precision x(n),fvec(n),diag(n),fjac(ldfjac,n),r(lr),
 
7
     *                 qtf(n),wa1(n),wa2(n),wa3(n),wa4(n)
 
8
      external fcn
 
9
c     **********
 
10
c
 
11
c     subroutine hybrd
 
12
c
 
13
c     the purpose of hybrd is to find a zero of a system of
 
14
c     n nonlinear functions in n variables by a modification
 
15
c     of the powell hybrid method. the user must provide a
 
16
c     subroutine which calculates the functions. the jacobian is
 
17
c     then calculated by a forward-difference approximation.
 
18
c
 
19
c     the subroutine statement is
 
20
c
 
21
c       subroutine hybrd(fcn,n,x,fvec,xtol,maxfev,ml,mu,epsfcn,
 
22
c                        diag,mode,factor,nprint,info,nfev,fjac,
 
23
c                        ldfjac,r,lr,qtf,wa1,wa2,wa3,wa4)
 
24
c
 
25
c     where
 
26
c
 
27
c       fcn is the name of the user-supplied subroutine which
 
28
c         calculates the functions. fcn must be declared
 
29
c         in an external statement in the user calling
 
30
c         program, and should be written as follows.
 
31
c
 
32
c         subroutine fcn(n,x,fvec,iflag)
 
33
c         integer n,iflag
 
34
c         double precision x(n),fvec(n)
 
35
c         ----------
 
36
c         calculate the functions at x and
 
37
c         return this vector in fvec.
 
38
c         ---------
 
39
c         return
 
40
c         end
 
41
c
 
42
c         the value of iflag should not be changed by fcn unless
 
43
c         the user wants to terminate execution of hybrd.
 
44
c         in this case set iflag to a negative integer.
 
45
c
 
46
c       n is a positive integer input variable set to the number
 
47
c         of functions and variables.
 
48
c
 
49
c       x is an array of length n. on input x must contain
 
50
c         an initial estimate of the solution vector. on output x
 
51
c         contains the final estimate of the solution vector.
 
52
c
 
53
c       fvec is an output array of length n which contains
 
54
c         the functions evaluated at the output x.
 
55
c
 
56
c       xtol is a nonnegative input variable. termination
 
57
c         occurs when the relative error between two consecutive
 
58
c         iterates is at most xtol.
 
59
c
 
60
c       maxfev is a positive integer input variable. termination
 
61
c         occurs when the number of calls to fcn is at least maxfev
 
62
c         by the end of an iteration.
 
63
c
 
64
c       ml is a nonnegative integer input variable which specifies
 
65
c         the number of subdiagonals within the band of the
 
66
c         jacobian matrix. if the jacobian is not banded, set
 
67
c         ml to at least n - 1.
 
68
c
 
69
c       mu is a nonnegative integer input variable which specifies
 
70
c         the number of superdiagonals within the band of the
 
71
c         jacobian matrix. if the jacobian is not banded, set
 
72
c         mu to at least n - 1.
 
73
c
 
74
c       epsfcn is an input variable used in determining a suitable
 
75
c         step length for the forward-difference approximation. this
 
76
c         approximation assumes that the relative errors in the
 
77
c         functions are of the order of epsfcn. if epsfcn is less
 
78
c         than the machine precision, it is assumed that the relative
 
79
c         errors in the functions are of the order of the machine
 
80
c         precision.
 
81
c
 
82
c       diag is an array of length n. if mode = 1 (see
 
83
c         below), diag is internally set. if mode = 2, diag
 
84
c         must contain positive entries that serve as
 
85
c         multiplicative scale factors for the variables.
 
86
c
 
87
c       mode is an integer input variable. if mode = 1, the
 
88
c         variables will be scaled internally. if mode = 2,
 
89
c         the scaling is specified by the input diag. other
 
90
c         values of mode are equivalent to mode = 1.
 
91
c
 
92
c       factor is a positive input variable used in determining the
 
93
c         initial step bound. this bound is set to the product of
 
94
c         factor and the euclidean norm of diag*x if nonzero, or else
 
95
c         to factor itself. in most cases factor should lie in the
 
96
c         interval (.1,100.). 100. is a generally recommended value.
 
97
c
 
98
c       nprint is an integer input variable that enables controlled
 
99
c         printing of iterates if it is positive. in this case,
 
100
c         fcn is called with iflag = 0 at the beginning of the first
 
101
c         iteration and every nprint iterations thereafter and
 
102
c         immediately prior to return, with x and fvec available
 
103
c         for printing. if nprint is not positive, no special calls
 
104
c         of fcn with iflag = 0 are made.
 
105
c
 
106
c       info is an integer output variable. if the user has
 
107
c         terminated execution, info is set to the (negative)
 
108
c         value of iflag. see description of fcn. otherwise,
 
109
c         info is set as follows.
 
110
c
 
111
c         info = 0   improper input parameters.
 
112
c
 
113
c         info = 1   relative error between two consecutive iterates
 
114
c                    is at most xtol.
 
115
c
 
116
c         info = 2   number of calls to fcn has reached or exceeded
 
117
c                    maxfev.
 
118
c
 
119
c         info = 3   xtol is too small. no further improvement in
 
120
c                    the approximate solution x is possible.
 
121
c
 
122
c         info = 4   iteration is not making good progress, as
 
123
c                    measured by the improvement from the last
 
124
c                    five jacobian evaluations.
 
125
c
 
126
c         info = 5   iteration is not making good progress, as
 
127
c                    measured by the improvement from the last
 
128
c                    ten iterations.
 
129
c
 
130
c       nfev is an integer output variable set to the number of
 
131
c         calls to fcn.
 
132
c
 
133
c       fjac is an output n by n array which contains the
 
134
c         orthogonal matrix q produced by the qr factorization
 
135
c         of the final approximate jacobian.
 
136
c
 
137
c       ldfjac is a positive integer input variable not less than n
 
138
c         which specifies the leading dimension of the array fjac.
 
139
c
 
140
c       r is an output array of length lr which contains the
 
141
c         upper triangular matrix produced by the qr factorization
 
142
c         of the final approximate jacobian, stored rowwise.
 
143
c
 
144
c       lr is a positive integer input variable not less than
 
145
c         (n*(n+1))/2.
 
146
c
 
147
c       qtf is an output array of length n which contains
 
148
c         the vector (q transpose)*fvec.
 
149
c
 
150
c       wa1, wa2, wa3, and wa4 are work arrays of length n.
 
151
c
 
152
c     subprograms called
 
153
c
 
154
c       user-supplied ...... fcn
 
155
c
 
156
c       minpack-supplied ... dogleg,dpmpar,enorm,fdjac1,
 
157
c                            qform,qrfac,r1mpyq,r1updt
 
158
c
 
159
c       fortran-supplied ... dabs,dmax1,dmin1,min0,mod
 
160
c
 
161
c     argonne national laboratory. minpack project. march 1980.
 
162
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
 
163
c
 
164
c     **********
 
165
      integer i,iflag,iter,j,jm1,l,msum,ncfail,ncsuc,nslow1,nslow2
 
166
      integer iwa(1)
 
167
      logical jeval,sing
 
168
      double precision actred,delta,epsmch,fnorm,fnorm1,one,pnorm,
 
169
     *                 prered,p1,p5,p001,p0001,ratio,sum,temp,xnorm,
 
170
     *                 zero
 
171
      double precision dpmpar,enorm
 
172
      data one,p1,p5,p001,p0001,zero
 
173
     *     /1.0d0,1.0d-1,5.0d-1,1.0d-3,1.0d-4,0.0d0/
 
174
c
 
175
c     epsmch is the machine precision.
 
176
c
 
177
      epsmch = dpmpar(1)
 
178
c
 
179
      info = 0
 
180
      iflag = 0
 
181
      nfev = 0
 
182
c
 
183
c     check the input parameters for errors.
 
184
c
 
185
      if (n .le. 0 .or. xtol .lt. zero .or. maxfev .le. 0
 
186
     *    .or. ml .lt. 0 .or. mu .lt. 0 .or. factor .le. zero
 
187
     *    .or. ldfjac .lt. n .or. lr .lt. (n*(n + 1))/2) go to 300
 
188
      if (mode .ne. 2) go to 20
 
189
      do 10 j = 1, n
 
190
         if (diag(j) .le. zero) go to 300
 
191
   10    continue
 
192
   20 continue
 
193
c
 
194
c     evaluate the function at the starting point
 
195
c     and calculate its norm.
 
196
c
 
197
      iflag = 1
 
198
      call fcn(n,x,fvec,iflag)
 
199
      nfev = 1
 
200
      if (iflag .lt. 0) go to 300
 
201
      fnorm = enorm(n,fvec)
 
202
c
 
203
c     determine the number of calls to fcn needed to compute
 
204
c     the jacobian matrix.
 
205
c
 
206
      msum = min0(ml+mu+1,n)
 
207
c
 
208
c     initialize iteration counter and monitors.
 
209
c
 
210
      iter = 1
 
211
      ncsuc = 0
 
212
      ncfail = 0
 
213
      nslow1 = 0
 
214
      nslow2 = 0
 
215
c
 
216
c     beginning of the outer loop.
 
217
c
 
218
   30 continue
 
219
         jeval = .true.
 
220
c
 
221
c        calculate the jacobian matrix.
 
222
c
 
223
         iflag = 2
 
224
         call fdjac1(fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn,wa1,
 
225
     *               wa2)
 
226
         nfev = nfev + msum
 
227
         if (iflag .lt. 0) go to 300
 
228
c
 
229
c        compute the qr factorization of the jacobian.
 
230
c
 
231
         call qrfac(n,n,fjac,ldfjac,.false.,iwa,1,wa1,wa2,wa3)
 
232
c
 
233
c        on the first iteration and if mode is 1, scale according
 
234
c        to the norms of the columns of the initial jacobian.
 
235
c
 
236
         if (iter .ne. 1) go to 70
 
237
         if (mode .eq. 2) go to 50
 
238
         do 40 j = 1, n
 
239
            diag(j) = wa2(j)
 
240
            if (wa2(j) .eq. zero) diag(j) = one
 
241
   40       continue
 
242
   50    continue
 
243
c
 
244
c        on the first iteration, calculate the norm of the scaled x
 
245
c        and initialize the step bound delta.
 
246
c
 
247
         do 60 j = 1, n
 
248
            wa3(j) = diag(j)*x(j)
 
249
   60       continue
 
250
         xnorm = enorm(n,wa3)
 
251
         delta = factor*xnorm
 
252
         if (delta .eq. zero) delta = factor
 
253
   70    continue
 
254
c
 
255
c        form (q transpose)*fvec and store in qtf.
 
256
c
 
257
         do 80 i = 1, n
 
258
            qtf(i) = fvec(i)
 
259
   80       continue
 
260
         do 120 j = 1, n
 
261
            if (fjac(j,j) .eq. zero) go to 110
 
262
            sum = zero
 
263
            do 90 i = j, n
 
264
               sum = sum + fjac(i,j)*qtf(i)
 
265
   90          continue
 
266
            temp = -sum/fjac(j,j)
 
267
            do 100 i = j, n
 
268
               qtf(i) = qtf(i) + fjac(i,j)*temp
 
269
  100          continue
 
270
  110       continue
 
271
  120       continue
 
272
c
 
273
c        copy the triangular factor of the qr factorization into r.
 
274
c
 
275
         sing = .false.
 
276
         do 150 j = 1, n
 
277
            l = j
 
278
            jm1 = j - 1
 
279
            if (jm1 .lt. 1) go to 140
 
280
            do 130 i = 1, jm1
 
281
               r(l) = fjac(i,j)
 
282
               l = l + n - i
 
283
  130          continue
 
284
  140       continue
 
285
            r(l) = wa1(j)
 
286
            if (wa1(j) .eq. zero) sing = .true.
 
287
  150       continue
 
288
c
 
289
c        accumulate the orthogonal factor in fjac.
 
290
c
 
291
         call qform(n,n,fjac,ldfjac,wa1)
 
292
c
 
293
c        rescale if necessary.
 
294
c
 
295
         if (mode .eq. 2) go to 170
 
296
         do 160 j = 1, n
 
297
            diag(j) = dmax1(diag(j),wa2(j))
 
298
  160       continue
 
299
  170    continue
 
300
c
 
301
c        beginning of the inner loop.
 
302
c
 
303
  180    continue
 
304
c
 
305
c           if requested, call fcn to enable printing of iterates.
 
306
c
 
307
            if (nprint .le. 0) go to 190
 
308
            iflag = 0
 
309
            if (mod(iter-1,nprint) .eq. 0) call fcn(n,x,fvec,iflag)
 
310
            if (iflag .lt. 0) go to 300
 
311
  190       continue
 
312
c
 
313
c           determine the direction p.
 
314
c
 
315
            call dogleg(n,r,lr,diag,qtf,delta,wa1,wa2,wa3)
 
316
c
 
317
c           store the direction p and x + p. calculate the norm of p.
 
318
c
 
319
            do 200 j = 1, n
 
320
               wa1(j) = -wa1(j)
 
321
               wa2(j) = x(j) + wa1(j)
 
322
               wa3(j) = diag(j)*wa1(j)
 
323
  200          continue
 
324
            pnorm = enorm(n,wa3)
 
325
c
 
326
c           on the first iteration, adjust the initial step bound.
 
327
c
 
328
            if (iter .eq. 1) delta = dmin1(delta,pnorm)
 
329
c
 
330
c           evaluate the function at x + p and calculate its norm.
 
331
c
 
332
            iflag = 1
 
333
            call fcn(n,wa2,wa4,iflag)
 
334
            nfev = nfev + 1
 
335
            if (iflag .lt. 0) go to 300
 
336
            fnorm1 = enorm(n,wa4)
 
337
c
 
338
c           compute the scaled actual reduction.
 
339
c
 
340
            actred = -one
 
341
            if (fnorm1 .lt. fnorm) actred = one - (fnorm1/fnorm)**2
 
342
c
 
343
c           compute the scaled predicted reduction.
 
344
c
 
345
            l = 1
 
346
            do 220 i = 1, n
 
347
               sum = zero
 
348
               do 210 j = i, n
 
349
                  sum = sum + r(l)*wa1(j)
 
350
                  l = l + 1
 
351
  210             continue
 
352
               wa3(i) = qtf(i) + sum
 
353
  220          continue
 
354
            temp = enorm(n,wa3)
 
355
            prered = zero
 
356
            if (temp .lt. fnorm) prered = one - (temp/fnorm)**2
 
357
c
 
358
c           compute the ratio of the actual to the predicted
 
359
c           reduction.
 
360
c
 
361
            ratio = zero
 
362
            if (prered .gt. zero) ratio = actred/prered
 
363
c
 
364
c           update the step bound.
 
365
c
 
366
            if (ratio .ge. p1) go to 230
 
367
               ncsuc = 0
 
368
               ncfail = ncfail + 1
 
369
               delta = p5*delta
 
370
               go to 240
 
371
  230       continue
 
372
               ncfail = 0
 
373
               ncsuc = ncsuc + 1
 
374
               if (ratio .ge. p5 .or. ncsuc .gt. 1)
 
375
     *            delta = dmax1(delta,pnorm/p5)
 
376
               if (dabs(ratio-one) .le. p1) delta = pnorm/p5
 
377
  240       continue
 
378
c
 
379
c           test for successful iteration.
 
380
c
 
381
            if (ratio .lt. p0001) go to 260
 
382
c
 
383
c           successful iteration. update x, fvec, and their norms.
 
384
c
 
385
            do 250 j = 1, n
 
386
               x(j) = wa2(j)
 
387
               wa2(j) = diag(j)*x(j)
 
388
               fvec(j) = wa4(j)
 
389
  250          continue
 
390
            xnorm = enorm(n,wa2)
 
391
            fnorm = fnorm1
 
392
            iter = iter + 1
 
393
  260       continue
 
394
c
 
395
c           determine the progress of the iteration.
 
396
c
 
397
            nslow1 = nslow1 + 1
 
398
            if (actred .ge. p001) nslow1 = 0
 
399
            if (jeval) nslow2 = nslow2 + 1
 
400
            if (actred .ge. p1) nslow2 = 0
 
401
c
 
402
c           test for convergence.
 
403
c
 
404
            if (delta .le. xtol*xnorm .or. fnorm .eq. zero) info = 1
 
405
            if (info .ne. 0) go to 300
 
406
c
 
407
c           tests for termination and stringent tolerances.
 
408
c
 
409
            if (nfev .ge. maxfev) info = 2
 
410
            if (p1*dmax1(p1*delta,pnorm) .le. epsmch*xnorm) info = 3
 
411
            if (nslow2 .eq. 5) info = 4
 
412
            if (nslow1 .eq. 10) info = 5
 
413
            if (info .ne. 0) go to 300
 
414
c
 
415
c           criterion for recalculating jacobian approximation
 
416
c           by forward differences.
 
417
c
 
418
            if (ncfail .eq. 2) go to 290
 
419
c
 
420
c           calculate the rank one modification to the jacobian
 
421
c           and update qtf if necessary.
 
422
c
 
423
            do 280 j = 1, n
 
424
               sum = zero
 
425
               do 270 i = 1, n
 
426
                  sum = sum + fjac(i,j)*wa4(i)
 
427
  270             continue
 
428
               wa2(j) = (sum - wa3(j))/pnorm
 
429
               wa1(j) = diag(j)*((diag(j)*wa1(j))/pnorm)
 
430
               if (ratio .ge. p0001) qtf(j) = sum
 
431
  280          continue
 
432
c
 
433
c           compute the qr factorization of the updated jacobian.
 
434
c
 
435
            call r1updt(n,n,r,lr,wa1,wa2,wa3,sing)
 
436
            call r1mpyq(n,n,fjac,ldfjac,wa2,wa3)
 
437
            call r1mpyq(1,n,qtf,1,wa2,wa3)
 
438
c
 
439
c           end of the inner loop.
 
440
c
 
441
            jeval = .false.
 
442
            go to 180
 
443
  290    continue
 
444
c
 
445
c        end of the outer loop.
 
446
c
 
447
         go to 30
 
448
  300 continue
 
449
c
 
450
c     termination, either normal or user imposed.
 
451
c
 
452
      if (iflag .lt. 0) info = iflag
 
453
      iflag = 0
 
454
      if (nprint .gt. 0) call fcn(n,x,fvec,iflag)
 
455
      return
 
456
c
 
457
c     last card of subroutine hybrd.
 
458
c
 
459
      end