1 |
heimbach |
1.2 |
|
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 |
heimbach |
1.3 |
external DDOT |
45 |
|
|
double precision DDOT |
46 |
heimbach |
1.2 |
|
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 |
heimbach |
1.3 |
r = DDOT( nn, dd, 1, xdiff,1 ) |
74 |
heimbach |
1.2 |
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 |
heimbach |
1.3 |
r = alpha(jp) - DDOT( nn, dd,1 , xdiff, 1 ) |
99 |
heimbach |
1.2 |
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 |