File: vecjac.f

package info (click to toggle)
cminpack 1.3.11-1
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 3,528 kB
  • sloc: ansic: 11,704; fortran: 5,648; makefile: 429; f90: 354; sh: 10
file content (321 lines) | stat: -rw-r--r-- 8,143 bytes parent folder | download | duplicates (7)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
      subroutine vecjac(n,x,fjac,ldfjac,nprob)
      integer n,ldfjac,nprob
      double precision x(n),fjac(ldfjac,n)
c     **********
c
c     subroutine vecjac
c
c     this subroutine defines the jacobian matrices of fourteen
c     test functions. the problem dimensions are as described
c     in the prologue comments of vecfcn.
c
c     the subroutine statement is
c
c       subroutine vecjac(n,x,fjac,ldfjac,nprob)
c
c     where
c
c       n is a positive integer variable.
c
c       x is an array of length n.
c
c       fjac is an n by n array. on output fjac contains the
c         jacobian matrix of the nprob function evaluated at x.
c
c       ldfjac is a positive integer variable not less than n
c         which specifies the leading dimension of the array fjac.
c
c       nprob is a positive integer variable which defines the
c         number of the problem. nprob must not exceed 14.
c
c     subprograms called
c
c       fortran-supplied ... datan,dcos,dexp,dmin1,dsin,dsqrt,
c                            max0,min0
c
c     argonne national laboratory. minpack project. march 1980.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     **********
      integer i,ivar,j,k,k1,k2,ml,mu
      double precision c1,c3,c4,c5,c6,c9,eight,fiftn,five,four,h,
     *                 hundrd,one,prod,six,sum,sum1,sum2,temp,temp1,
     *                 temp2,temp3,temp4,ten,three,ti,tj,tk,tpi,
     *                 twenty,two,zero
      double precision dfloat
      data zero,one,two,three,four,five,six,eight,ten,fiftn,twenty,
     *     hundrd
     *     /0.0d0,1.0d0,2.0d0,3.0d0,4.0d0,5.0d0,6.0d0,8.0d0,1.0d1,
     *      1.5d1,2.0d1,1.0d2/
      data c1,c3,c4,c5,c6,c9 /1.0d4,2.0d2,2.02d1,1.98d1,1.8d2,2.9d1/
      dfloat(ivar) = ivar
c
c     jacobian routine selector.
c
      go to (10,20,50,60,90,100,200,230,290,320,350,380,420,450),
     *      nprob
c
c     rosenbrock function.
c
   10 continue
      fjac(1,1) = -one
      fjac(1,2) = zero
      fjac(2,1) = -twenty*x(1)
      fjac(2,2) = ten
      go to 490
c
c     powell singular function.
c
   20 continue
      do 40 k = 1, 4
         do 30 j = 1, 4
            fjac(k,j) = zero
   30       continue
   40    continue
      fjac(1,1) = one
      fjac(1,2) = ten
      fjac(2,3) = dsqrt(five)
      fjac(2,4) = -fjac(2,3)
      fjac(3,2) = two*(x(2) - two*x(3))
      fjac(3,3) = -two*fjac(3,2)
      fjac(4,1) = two*dsqrt(ten)*(x(1) - x(4))
      fjac(4,4) = -fjac(4,1)
      go to 490
c
c     powell badly scaled function.
c
   50 continue
      fjac(1,1) = c1*x(2)
      fjac(1,2) = c1*x(1)
      fjac(2,1) = -dexp(-x(1))
      fjac(2,2) = -dexp(-x(2))
      go to 490
c
c     wood function.
c
   60 continue
      do 80 k = 1, 4
         do 70 j = 1, 4
            fjac(k,j) = zero
   70       continue
   80    continue
      temp1 = x(2) - three*x(1)**2
      temp2 = x(4) - three*x(3)**2
      fjac(1,1) = -c3*temp1 + one
      fjac(1,2) = -c3*x(1)
      fjac(2,1) = -two*c3*x(1)
      fjac(2,2) = c3 + c4
      fjac(2,4) = c5
      fjac(3,3) = -c6*temp2 + one
      fjac(3,4) = -c6*x(3)
      fjac(4,2) = c5
      fjac(4,3) = -two*c6*x(3)
      fjac(4,4) = c6 + c4
      go to 490
c
c     helical valley function.
c
   90 continue
      tpi = eight*datan(one)
      temp = x(1)**2 + x(2)**2
      temp1 = tpi*temp
      temp2 = dsqrt(temp)
      fjac(1,1) = hundrd*x(2)/temp1
      fjac(1,2) = -hundrd*x(1)/temp1
      fjac(1,3) = ten
      fjac(2,1) = ten*x(1)/temp2
      fjac(2,2) = ten*x(2)/temp2
      fjac(2,3) = zero
      fjac(3,1) = zero
      fjac(3,2) = zero
      fjac(3,3) = one
      go to 490
c
c     watson function.
c
  100 continue
      do 120 k = 1, n
         do 110 j = k, n
            fjac(k,j) = zero
  110       continue
  120    continue
      do 170 i = 1, 29
         ti = dfloat(i)/c9
         sum1 = zero
         temp = one
         do 130 j = 2, n
            sum1 = sum1 + dfloat(j-1)*temp*x(j)
            temp = ti*temp
  130       continue
         sum2 = zero
         temp = one
         do 140 j = 1, n
            sum2 = sum2 + temp*x(j)
            temp = ti*temp
  140       continue
         temp1 = two*(sum1 - sum2**2 - one)
         temp2 = two*sum2
         temp = ti**2
         tk = one
         do 160 k = 1, n
            tj = tk
            do 150 j = k, n
               fjac(k,j) = fjac(k,j)
     *                     + tj
     *                       *((dfloat(k-1)/ti - temp2)
     *                         *(dfloat(j-1)/ti - temp2) - temp1)
               tj = ti*tj
  150          continue
            tk = temp*tk
  160       continue
  170    continue
      fjac(1,1) = fjac(1,1) + six*x(1)**2 - two*x(2) + three
      fjac(1,2) = fjac(1,2) - two*x(1)
      fjac(2,2) = fjac(2,2) + one
      do 190 k = 1, n
         do 180 j = k, n
            fjac(j,k) = fjac(k,j)
  180       continue
  190    continue
      go to 490
c
c     chebyquad function.
c
  200 continue
      tk = one/dfloat(n)
      do 220 j = 1, n
         temp1 = one
         temp2 = two*x(j) - one
         temp = two*temp2
         temp3 = zero
         temp4 = two
         do 210 k = 1, n
            fjac(k,j) = tk*temp4
            ti = four*temp2 + temp*temp4 - temp3
            temp3 = temp4
            temp4 = ti
            ti = temp*temp2 - temp1
            temp1 = temp2
            temp2 = ti
  210       continue
  220    continue
      go to 490
c
c     brown almost-linear function.
c
  230 continue
      prod = one
      do 250 j = 1, n
         prod = x(j)*prod
         do 240 k = 1, n
            fjac(k,j) = one
  240       continue
         fjac(j,j) = two
  250    continue
      do 280 j = 1, n
         temp = x(j)
         if (temp .ne. zero) go to 270
         temp = one
         prod = one
         do 260 k = 1, n
            if (k .ne. j) prod = x(k)*prod
  260       continue
  270    continue
         fjac(n,j) = prod/temp
  280    continue
      go to 490
c
c     discrete boundary value function.
c
  290 continue
      h = one/dfloat(n+1)
      do 310 k = 1, n
         temp = three*(x(k) + dfloat(k)*h + one)**2
         do 300 j = 1, n
            fjac(k,j) = zero
  300       continue
         fjac(k,k) = two + temp*h**2/two
         if (k .ne. 1) fjac(k,k-1) = -one
         if (k .ne. n) fjac(k,k+1) = -one
  310    continue
      go to 490
c
c     discrete integral equation function.
c
  320 continue
      h = one/dfloat(n+1)
      do 340 k = 1, n
         tk = dfloat(k)*h
         do 330 j = 1, n
            tj = dfloat(j)*h
            temp = three*(x(j) + tj + one)**2
            fjac(k,j) = h*dmin1(tj*(one-tk),tk*(one-tj))*temp/two
  330       continue
         fjac(k,k) = fjac(k,k) + one
  340    continue
      go to 490
c
c     trigonometric function.
c
  350 continue
      do 370 j = 1, n
         temp = dsin(x(j))
         do 360 k = 1, n
            fjac(k,j) = temp
  360       continue
         fjac(j,j) = dfloat(j+1)*temp - dcos(x(j))
  370    continue
      go to 490
c
c     variably dimensioned function.
c
  380 continue
      sum = zero
      do 390 j = 1, n
         sum = sum + dfloat(j)*(x(j) - one)
  390    continue
      temp = one + six*sum**2
      do 410 k = 1, n
         do 400 j = k, n
            fjac(k,j) = dfloat(k*j)*temp
            fjac(j,k) = fjac(k,j)
  400       continue
         fjac(k,k) = fjac(k,k) + one
  410    continue
      go to 490
c
c     broyden tridiagonal function.
c
  420 continue
      do 440 k = 1, n
         do 430 j = 1, n
            fjac(k,j) = zero
  430       continue
         fjac(k,k) = three - four*x(k)
         if (k .ne. 1) fjac(k,k-1) = -one
         if (k .ne. n) fjac(k,k+1) = -two
  440    continue
      go to 490
c
c     broyden banded function.
c
  450 continue
      ml = 5
      mu = 1
      do 480 k = 1, n
         do 460 j = 1, n
            fjac(k,j) = zero
  460       continue
         k1 = max0(1,k-ml)
         k2 = min0(k+mu,n)
         do 470 j = k1, k2
            if (j .ne. k) fjac(k,j) = -(one + two*x(j))
  470       continue
         fjac(k,k) = two + fiftn*x(k)**2
  480    continue
  490 continue
      return
c
c     last card of subroutine vecjac.
c
      end