/[ascend]/trunk/minos54/mi65rmod.f
ViewVC logotype

Contents of /trunk/minos54/mi65rmod.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (show annotations) (download)
Fri Oct 29 20:54:12 2004 UTC (19 years, 7 months ago) by aw0a
File size: 16666 byte(s)
Setting up web subdirectory in repository
1 ************************************************************************
2 *
3 * File mi65rmod fortran.
4 *
5 * m6bfgs m6bswp m6radd m6rcnd m6rdel
6 * m6rmod m6rset m6rsol m6swap
7 *
8 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
9
10 subroutine m6bfgs( maxr, n, nr, r, g, g2, p, v,
11 $ step, told, tolz, inform )
12
13 implicit double precision (a-h,o-z)
14 double precision r(nr), g(n), g2(n), p(n), v(n)
15
16 * ------------------------------------------------------------------
17 * m6bfgs applies the BFGS update to the upper-triangular matrix r,
18 * which holds the cholesky factor of the quasi-newton approximation
19 * of the (projected) hessian.
20 *
21 * r contains a triangle of size ncolr = min( n, maxr ).
22 * If n .gt. maxr, r also contains a diagonal of size n - maxr.
23 *
24 * p holds the search direction. It is overwritten.
25 * v must satisfy r(t)*v = g. It is overwritten.
26 *
27 * On exit,
28 * inform = 0 if no update was performed,
29 * = 1 if the update was successful,
30 * = 2 if it was nearly singular.
31 * ------------------------------------------------------------------
32
33 intrinsic abs, min, sqrt
34 parameter ( one = 1.0d+0 )
35
36 inform = 0
37 ncolr = min( n, maxr )
38 gtp = ddot ( n, g , 1, p, 1 )
39 gtp2 = ddot ( n, g2, 1, p, 1 )
40 if (gtp2 .le. 0.91*gtp) return
41
42 delta1 = one / sqrt( abs( gtp ) )
43 delta2 = one / sqrt( step*(gtp2 - gtp) )
44
45 * Normalize v and change its sign.
46
47 call dscal ( n, (- delta1), v, 1 )
48
49 * Avoid cancellation error in forming the new vector p.
50
51 if ( abs( delta1/delta2 - one ) .ge. 0.5d+0) then
52 do 100 j = 1, n
53 p(j) = delta2*( g2(j) - g(j) ) + delta1*g(j)
54 100 continue
55 else
56 d = delta1 - delta2
57 do 120 j = 1, n
58 p(j) = delta2*g2(j) + d*g(j)
59 120 continue
60 end if
61
62 * Find the last nonzero in v.
63
64 do 450 ii = 1, ncolr
65 lastv = ncolr + 1 - ii
66 if (abs( v(lastv) ) .gt. tolz) go to 500
67 450 continue
68 lastv = 1
69
70 * Triangularize r + v * p(t).
71
72 500 call m6rmod( ncolr, nr, r, v, p, lastv, told, tolz, inform )
73
74 * Deal with surplus diagonals of r.
75
76 if (n .gt. maxr) then
77 j1 = maxr + 1
78 l = maxr*j1/2
79 do 650 j = j1, n
80 l = l + 1
81 r(l) = r(l) + v(j)*p(j)
82 650 continue
83 end if
84
85 * end of m6bfgs
86 end
87
88 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
89
90 subroutine m6bswp( n, nr, r, v, w, lastv, told, tolz, inform )
91
92 implicit double precision (a-h,o-z)
93 double precision r(nr), v(n), w(n)
94
95 * ------------------------------------------------------------------
96 * m6bswp modifies the upper-triangular matrix r
97 * to account for a basis exchange in which the lastv-th
98 * superbasic variable becomes basic. r is changed to
99 * r + v*w(t), which is triangularized by r1mod,
100 * where v is the lastv-th column of r, and w is input.
101 * ------------------------------------------------------------------
102
103 * Set v = lastv-th column of r and find its norm.
104
105 l = (lastv - 1)*lastv/2
106 call dcopy ( lastv, r(l+1), 1, v, 1 )
107 vnorm = dasum ( lastv, v, 1 )
108 call m6rmod( n, nr, r, v, w, lastv,
109 $ (vnorm*told), (vnorm*tolz), inform )
110
111 * end of m6bswp
112 end
113
114 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
115
116 subroutine m6radd( maxr, nr, ns, r )
117
118 implicit double precision (a-h,o-z)
119 double precision r(nr)
120
121 * ------------------------------------------------------------------
122 * m6radd adds column ns to the upper triangular matrix r.
123 * In this version (Sep 1984) it is just a unit vector.
124 *
125 * Modified Jan 1983 to add a diagonal only, if ns .gt. maxr.
126 * ------------------------------------------------------------------
127
128 parameter ( zero = 0.0d+0, one = 1.0d+0 )
129
130 if (ns .eq. 1) then
131 r(1) = one
132 else if (ns .le. maxr) then
133 ldiag1 = (ns - 1)*ns/2
134 ldiag2 = ldiag1 + ns
135 call dload ( ns, zero, r(ldiag1+1), 1 )
136 r(ldiag2) = one
137 else
138 ldiag2 = maxr*(maxr + 1)/2 + (ns - maxr)
139 r(ldiag2) = one
140 end if
141
142 * end of m6radd
143 end
144
145 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
146
147 subroutine m6rcnd( maxr, nr, ns, r, dmax, dmin, cond )
148
149 implicit double precision (a-h,o-z)
150 double precision r(nr)
151
152 * ------------------------------------------------------------------
153 * m6rcnd finds the largest and smallest diagonals of the
154 * upper triangular matrix r, and returns the square of their ratio.
155 * This is a lower bound on the condition number of r(t)*r.
156 * ------------------------------------------------------------------
157
158 intrinsic abs, max, min
159
160 ncolr = min( maxr, ns )
161 dmax = abs( r(1) )
162 dmin = dmax
163
164 if (ncolr .gt. 1) then
165 l = 1
166 do 100 j = 2, ncolr
167 l = l + j
168 d = abs( r(l) )
169 dmax = max( dmax, d )
170 dmin = min( dmin, d )
171 100 continue
172 end if
173
174 cond = (dmax/dmin)**2
175
176 * end of m6rcnd
177 end
178
179 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
180
181 subroutine m6rdel( m, maxr, nr, ns, ms,
182 $ kb, bbl, bbu, grd, r, rg, x, jq, rset )
183
184 implicit double precision (a-h,o-z)
185 integer kb(ms)
186 double precision bbl(ms), bbu(ms), grd(ms)
187 double precision r(nr), rg(ns), x(ms)
188 logical rset
189
190 * ------------------------------------------------------------------
191 * m6rdel deletes the jq-th superbasic variable from r
192 * and from various arrays kb, bbl, bbu, grd, rg, x.
193 * ------------------------------------------------------------------
194
195 common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy
196
197 intrinsic abs, max, min, sqrt
198 parameter ( zero = 0.0d+0 )
199
200 if (jq .eq. ns) return
201 if (.not. rset) go to 500
202
203 * Delete the jq-th column of r and triangularize the
204 * remaining columns using a partial forward sweep of rotations.
205
206 ncolr = min( maxr, ns )
207 k = jq*(jq + 1)/2
208
209 do 50 i = jq + 1, ncolr
210 k = k + i
211 t1 = r(k-1)
212 t2 = r(k)
213 d = sqrt( t1*t1 + t2*t2 )
214 r(k-1) = d
215 if (i .eq. ncolr) go to 20
216 sn = t2/d
217 cs = t1/d
218 if (abs( cs ) .le. eps0) cs = zero
219 j1 = i + 1
220 k1 = k + i
221 do 10 j = j1, ncolr
222 t1 = r(k1-1)
223 t2 = r(k1)
224 r(k1-1) = cs*t1 + sn*t2
225 r(k1) = sn*t1 - cs*t2
226 k1 = k1 + j
227 10 continue
228 20 k1 = i - 1
229 j2 = k - i
230 j1 = j2 - i + 2
231 do 30 j = j1, j2
232 r(j) = r(j+k1)
233 30 continue
234 50 continue
235
236 * If necessary, clear out the last column of r.
237
238 if (ns .gt. maxr) then
239 lastr = maxr*(maxr + 1)/2
240 l = lastr - maxr + 1
241 if (jq .le. maxr) call dload ( maxr, zero, r(l), 1 )
242
243 * Shift surplus diagonals of r to the left.
244
245 j1 = max( maxr, jq )
246 j2 = ns - 1
247 l = lastr + (j1 - maxr)
248 do 320 j = j1, j2
249 r(l) = r(l+1)
250 320 continue
251 end if
252
253 * Shift all the arrays one place to the left.
254
255 500 j2 = ns - 1
256 do 550 j = jq, j2
257 rg(j) = rg(j+1)
258 k = m + j
259 kb(k) = kb(k+1)
260 bbl(k) = bbl(k+1)
261 bbu(k) = bbu(k+1)
262 grd(k) = grd(k+1)
263 x(k) = x(k+1)
264 550 continue
265
266 * end of m6rdel
267 end
268
269 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
270
271 subroutine m6rmod( n, nr, r, v, w, lastv, told, tolz, inform )
272
273 implicit double precision (a-h,o-z)
274 double precision r(nr), v(n), w(n)
275
276 * ------------------------------------------------------------------
277 * m6rmod modifies the upper-triangular matrix r so that
278 * q*(r + v*w(t)) is upper triangular, where q is orthogonal,
279 * v and w are vectors, and the new r overwrites the old.
280 *
281 * q is the product of two sweeps of plane rotations (not stored).
282 * These affect the lastv-th row of r, which is temporarily held
283 * in the vector v. Thus, v is overwritten. w is not altered.
284 *
285 * lastv points to the last nonzero of v. The value lastv = n
286 * would always be ok, but sometimes it is known to be less
287 * than n, so q reduces to two partial sweeps of
288 * rotations.
289 *
290 * told is a tolerance on the lastv-th diagonal of r.
291 * tolz is a tolerance for negligible elements in v.
292 *
293 * On exit,
294 * inform = 1 if the diagonal of r is larger than told,
295 * = 2 if not (in which case it is reset to told).
296 * ------------------------------------------------------------------
297
298 intrinsic sqrt
299 parameter ( zero = 0.0d+0 )
300
301 vlast = v(lastv)
302 lm1 = lastv - 1
303 lastr = lastv*(lastv + 1)/2
304
305 * Copy the lastv-th row of r into the end of v.
306
307 l = lastr
308 do 100 j = lastv, n
309 v(j) = r(l)
310 l = l + j
311 100 continue
312
313 * Reduce v to a multiple of e(lastv)
314 * using a partial backward sweep of rotations.
315 * This fills in the lastv-th row of r (held in v).
316
317 v2 = vlast**2
318 ldiag = lastr
319
320 do 300 ii = 1, lm1
321 i = lastv - ii
322 ldiag = ldiag - i - 1
323 s = v(i)
324 v(i) = zero
325 if (abs(s) .le. tolz) go to 300
326 v2 = s**2 + v2
327 root = sqrt(v2)
328 cs = vlast/root
329 sn = s/root
330 vlast = root
331 l = ldiag
332
333 do 200 j = i, n
334 s = v(j)
335 t = r(l)
336 v(j) = cs*s + sn*t
337 r(l) = sn*s - cs*t
338 l = l + j
339 200 continue
340 300 continue
341
342 * Set v = v + vlast*w.
343
344 call daxpy ( n, vlast, w, 1, v, 1 )
345
346 * Eliminate the front of the lastv-th row of r (held in v)
347 * using a partial forward sweep of rotations.
348
349 ldiag = 0
350
351 do 700 i = 1, lm1
352 ldiag = ldiag + i
353 t = v(i)
354 if (abs(t) .le. tolz) go to 700
355 s = r(ldiag)
356 root = sqrt(s**2 + t**2)
357 cs = s/root
358 sn = t/root
359 r(ldiag) = root
360 l = ldiag + i
361
362 do 600 j = i + 1, n
363 s = r(l)
364 t = v(j)
365 r(l) = cs*s + sn*t
366 v(j) = sn*s - cs*t
367 l = l + j
368 600 continue
369 700 continue
370
371 * Insert the new lastv-th row of r.
372
373 l = lastr
374 do 900 j = lastv, n
375 r(l) = v(j)
376 l = l + j
377 900 continue
378
379 * Test for (unlikely) singularity.
380
381 inform = 1
382 if (abs( r(lastr) ) .le. told) then
383 inform = 2
384 r(lastr) = told
385 end if
386
387 * end of m6rmod
388 end
389
390 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
391
392 subroutine m6rset( maxr, nr, ns, r, cond )
393
394 implicit double precision (a-h,o-z)
395 double precision r(nr)
396
397 * ------------------------------------------------------------------
398 * m6rset alters r, the upper-triangular factor of the
399 * approximation to the reduced hessian.
400 *
401 * If r(1) = zero, r does not exist.
402 * In this case, r is initialized to the identity matrix.
403 *
404 * Otherwise, r already exists and we attempt to make it better
405 * conditioned by scaling its columns by the square roots of its
406 * current diagonals.
407 * ------------------------------------------------------------------
408
409 common /m1file/ iread,iprint,isumm
410
411 intrinsic abs, max, min, sqrt
412 parameter ( zero = 0.0d+0, one = 1.0d+0 )
413
414 cond = one
415 ncolr = min( maxr, ns )
416 if (ncolr .eq. 0 ) return
417
418 dmax = abs( r(1) )
419 if (dmax .eq. zero) then
420
421 * Set r = the identity.
422
423 l = ncolr*(ncolr + 1)/2
424 call dload ( l, zero, r, 1 )
425 ldiag = 0
426
427 do 100 k = 1, ncolr
428 ldiag = ldiag + k
429 r(ldiag) = one
430 100 continue
431
432 do 120 k = maxr + 1, ns
433 ldiag = ldiag + 1
434 r(ldiag) = one
435 120 continue
436 else
437
438 * Scale the columns of r.
439
440 dmin = dmax
441 ldiag = 0
442
443 do 250 k = 1, ncolr
444 l = ldiag + 1
445 ldiag = ldiag + k
446 diag = abs( r(ldiag) )
447 dmax = max( dmax, diag )
448 dmin = min( dmin, diag )
449 s = one / sqrt( diag )
450 do 240 i = l, ldiag
451 r(i) = r(i)*s
452 240 continue
453 250 continue
454
455 cond = dmax / dmin
456 if (iprint .gt. 0) write(iprint, 2000) dmax, dmin, cond
457 end if
458 return
459
460 2000 format(/ ' XXX Hessian modified. Diags =', 1p, 2e11.2,
461 $ ', Cond R(t)*R =', e11.2)
462
463 * end of m6rset
464 end
465
466 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
467
468 subroutine m6rsol( mode, maxr, nr, ns, r, y )
469
470 implicit double precision (a-h,o-z)
471 double precision r(nr), y(ns)
472
473 * ------------------------------------------------------------------
474 * m6rsol solves R*x = y or R(t)*x = y, where R is an
475 * upper-triangular matrix stored by columns in r(nr).
476 * The solution x overwrites y.
477 * ------------------------------------------------------------------
478
479 intrinsic min
480
481 ncolr = min( maxr, ns )
482
483 if (mode .eq. 1) then
484
485 * mode = 1 -- solve R*y = y.
486
487 l = ncolr*(ncolr + 1)/2
488 do 120 i = ncolr, 2, -1
489 t = y(i) / r(l)
490 y(i) = t
491 l = l - i
492 call daxpy ( i-1, (-t), r(l+1), 1, y, 1 )
493 120 continue
494 y(1) = y(1) / r(1)
495 else
496
497 * mode = 2 -- solve R(t)*y = y.
498
499 y(1) = y(1) / r(1)
500 l = 1
501 do 220 i = 2, ncolr
502 t = ddot ( i-1, r(l+1), 1, y, 1 )
503 l = l + i
504 y(i) = (y(i) - t) / r(l)
505 220 continue
506 end if
507
508 * Deal with surplus diagonals of r.
509
510 if (ns .gt. maxr) then
511 l = maxr*i/2
512 do 320 j = maxr + 1, ns
513 l = l + 1
514 y(j) = y(j) / r(l)
515 320 continue
516 end if
517
518 * end of m6rsol
519 end
520
521 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
522
523 subroutine m6swap( m, maxr, nr, ns, ms,
524 $ kb, bbl, bbu, grd, r, rg, x )
525
526 implicit double precision (a-h,o-z)
527 integer kb(ms)
528 double precision bbl(ms), bbu(ms), grd(ms)
529 double precision r(nr), rg(ns), x(ms)
530
531 * ------------------------------------------------------------------
532 * m6swap (superbasic swap) finds the largest reduced gradient in
533 * the range rg(maxr+1), ..., rg(ns) and swaps it into position
534 * maxr + 1 (so we know where it is).
535 * ------------------------------------------------------------------
536
537 k1 = maxr + 1
538 if (ns .gt. k1) then
539 nz2 = ns - maxr
540 k2 = maxr + idamax( nz2, rg(k1), 1 )
541 if (k2 .gt. k1) then
542 j = m + k1
543 k = m + k2
544 lastr = maxr*k1/2
545 ldiag1 = lastr + 1
546 ldiag2 = lastr + (k2 - maxr)
547
548 rdiag1 = r(ldiag1)
549 rg1 = rg(k1)
550 j1 = kb(j)
551 bl1 = bbl(j)
552 bu1 = bbu(j)
553 grd1 = grd(j)
554 x1 = x(j)
555
556 r(ldiag1) = r(ldiag2)
557 rg(k1) = rg(k2)
558 kb(j) = kb(k)
559 bbl(j) = bbl(k)
560 bbu(j) = bbu(k)
561 grd(j) = grd(k)
562 x(j) = x(k)
563
564 r(ldiag2) = rdiag1
565 rg(k2) = rg1
566 kb(k) = j1
567 bbl(k) = bl1
568 bbu(k) = bu1
569 grd(k) = grd1
570 x(k) = x1
571 end if
572 end if
573
574 * end of m6swap
575 end

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22