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

Annotation of /trunk/minos54/mi65rmod.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (hide 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 aw0a 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