1 |
|
2 |
subroutine hessupd( nn, mupd, dd, jmin, jmax, xdiff, lphprint ) |
3 |
|
4 |
c ================================================================== |
5 |
c SUBROUTINE hessupd |
6 |
c ================================================================== |
7 |
c |
8 |
c o controls update of descent vector using available |
9 |
c approximation of Hessian Matrix based on gradients of |
10 |
c previous iterations |
11 |
c |
12 |
c o Reference: J.C. Gilbert & C. Lemarechal |
13 |
c Some numerical experiments with variable-storage |
14 |
c quasi-Newton algorithms |
15 |
c Mathematical Programming 45 (1989), pp. 407-435 |
16 |
c |
17 |
c o started: ??? not reproducible |
18 |
c |
19 |
c o changed: Patrick Heimbach, MIT/EAPS |
20 |
c 24-Feb-2000: |
21 |
c - changed some variable names to be consistent |
22 |
c with routine lsoptv, lsline; |
23 |
c |
24 |
c o Version: 2.1.0, 02-Mar-2000: Patrick Heimbach, MIT/EAPS |
25 |
c |
26 |
c ================================================================== |
27 |
c SUBROUTINE hessupd |
28 |
c ================================================================== |
29 |
|
30 |
implicit none |
31 |
|
32 |
#include "blas1.h" |
33 |
|
34 |
c------------------------------------ |
35 |
c declare arguments |
36 |
c------------------------------------ |
37 |
integer nn, mupd, jmin, jmax |
38 |
double precision dd(nn), alpha(100), xdiff(nn) |
39 |
logical lphprint |
40 |
|
41 |
c------------------------------------ |
42 |
c declare local variables |
43 |
c------------------------------------ |
44 |
external DDOT |
45 |
double precision DDOT |
46 |
|
47 |
integer jfin, i, j, jp |
48 |
double precision r |
49 |
|
50 |
c------------------------------------ |
51 |
c initialization |
52 |
c------------------------------------ |
53 |
jfin = jmax |
54 |
|
55 |
if (lphprint) |
56 |
& print *, 'pathei-lsopt: in hessupd; ', |
57 |
& 'jmin, jmax, mupd:', jmin, jmax, mupd |
58 |
|
59 |
if (jfin.lt.jmin) jfin = jmax+mupd |
60 |
|
61 |
c------------------------------------ |
62 |
c compute right hand side |
63 |
c------------------------------------ |
64 |
do j = jfin,jmin,-1 |
65 |
|
66 |
if (lphprint) |
67 |
& print *, 'pathei-lsopt: in hessupd; loop ', |
68 |
& 'j,jfin,jmin = ', j,jfin,jmin |
69 |
|
70 |
jp = j |
71 |
if (jp.gt.mupd) jp = jp-mupd |
72 |
C xdiff is control difference (x(i+1)-x(i)) from previous iterations. |
73 |
call dostore( nn, xdiff, .false., 2*jp+3 ) |
74 |
r = DDOT( nn, dd, 1, xdiff,1 ) |
75 |
C xdiff now is gradient difference from previous iterations. |
76 |
call dostore( nn, xdiff, .false., 2*jp+2 ) |
77 |
alpha(jp) = r |
78 |
do i = 1, nn |
79 |
dd(i) = dd(i) - r*xdiff(i) |
80 |
end do |
81 |
end do |
82 |
|
83 |
c------------------------------------ |
84 |
c multiply precondition matrix |
85 |
c------------------------------------ |
86 |
if (mupd .ne. 0) then |
87 |
call dostore( nn, xdiff, .false., 3 ) |
88 |
do i = 1, nn |
89 |
dd(i) = dd(i)*xdiff(i) |
90 |
end do |
91 |
end if |
92 |
|
93 |
c------------------------------------ |
94 |
c compute left hand side |
95 |
c------------------------------------ |
96 |
do j = jmin,jfin |
97 |
jp = j |
98 |
if (jp .gt. mupd) jp = jp-mupd |
99 |
call dostore( nn, xdiff, .false., 2*jp+2 ) |
100 |
r = alpha(jp) - DDOT( nn, dd,1 , xdiff, 1 ) |
101 |
call dostore( nn, xdiff, .false., 2*jp+3 ) |
102 |
do i = 1, nn |
103 |
dd(i) = dd(i) + r*xdiff(i) |
104 |
end do |
105 |
end do |
106 |
|
107 |
return |
108 |
|
109 |
end |