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 |
call dostore( nn, xdiff, .false., 2*jp+3 ) |
73 |
r = DDOT( nn, dd, 1, xdiff,1 ) |
74 |
call dostore( nn, xdiff, .false., 2*jp+2 ) |
75 |
alpha(jp) = r |
76 |
do i = 1, nn |
77 |
dd(i) = dd(i) - r*xdiff(i) |
78 |
end do |
79 |
end do |
80 |
|
81 |
c------------------------------------ |
82 |
c multiply precondition matrix |
83 |
c------------------------------------ |
84 |
if (mupd .ne. 0) then |
85 |
call dostore( nn, xdiff, .false., 3 ) |
86 |
do i = 1, nn |
87 |
dd(i) = dd(i)*xdiff(i) |
88 |
end do |
89 |
end if |
90 |
|
91 |
c------------------------------------ |
92 |
c compute left hand side |
93 |
c------------------------------------ |
94 |
do j = jmin,jfin |
95 |
jp = j |
96 |
if (jp .gt. mupd) jp = jp-mupd |
97 |
call dostore( nn, xdiff, .false., 2*jp+2 ) |
98 |
r = alpha(jp) - DDOT( nn, dd,1 , xdiff, 1 ) |
99 |
call dostore( nn, xdiff, .false., 2*jp+3 ) |
100 |
do i = 1, nn |
101 |
dd(i) = dd(i) + r*xdiff(i) |
102 |
end do |
103 |
end do |
104 |
|
105 |
return |
106 |
|
107 |
end |