1 |
heimbach |
1.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 |
|
|
include 'mpif.h' |
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 |
|
|
integer status(MPI_STATUS_SIZE),ierr |
41 |
|
|
integer myid, nprocs,mystart,myend |
42 |
|
|
common /mpi_parm/ myid,nprocs,mystart,myend |
43 |
|
|
|
44 |
|
|
c------------------------------------ |
45 |
|
|
c declare local variables |
46 |
|
|
c------------------------------------ |
47 |
|
|
external DDOT |
48 |
|
|
double precision DDOT |
49 |
|
|
|
50 |
|
|
integer jfin, i, j, jp |
51 |
|
|
double precision r |
52 |
|
|
double precision temp |
53 |
|
|
|
54 |
|
|
c------------------------------------ |
55 |
|
|
c initialization |
56 |
|
|
c------------------------------------ |
57 |
|
|
jfin = jmax |
58 |
|
|
|
59 |
|
|
if (lphprint .and. myid .eq. 0) |
60 |
|
|
& print *, 'pathei-lsopt: in hessupd; ', |
61 |
|
|
& 'jmin, jmax, mupd:', jmin, jmax, mupd |
62 |
|
|
|
63 |
|
|
|
64 |
|
|
if (jfin.lt.jmin) jfin = jmax+mupd |
65 |
|
|
|
66 |
|
|
c------------------------------------ |
67 |
|
|
c compute right hand side |
68 |
|
|
c------------------------------------ |
69 |
|
|
do j = jfin,jmin,-1 |
70 |
|
|
|
71 |
|
|
if (lphprint .and. myid .eq. 0) |
72 |
|
|
& print *, 'pathei-lsopt: in hessupd; loop ', |
73 |
|
|
& 'j,jfin,jmin = ', j,jfin,jmin |
74 |
|
|
|
75 |
|
|
jp = j |
76 |
|
|
if (jp.gt.mupd) jp = jp-mupd |
77 |
|
|
call dostore( nn, xdiff, .false., 2*jp+3 ) |
78 |
|
|
|
79 |
|
|
temp = DDOT( nn, dd, 1, xdiff, 1 ) |
80 |
|
|
call MPI_ALLREDUCE(temp,r,1,MPI_DOUBLE_PRECISION,MPI_SUM, |
81 |
|
|
& MPI_COMM_WORLD,ierr) |
82 |
|
|
! r = DDOT( nn, dd, 1, xdiff,1 ) |
83 |
|
|
|
84 |
|
|
call dostore( nn, xdiff, .false., 2*jp+2 ) |
85 |
|
|
alpha(jp) = r |
86 |
|
|
do i = 1, nn |
87 |
|
|
dd(i) = dd(i) - r*xdiff(i) |
88 |
|
|
end do |
89 |
|
|
end do |
90 |
|
|
|
91 |
|
|
c------------------------------------ |
92 |
|
|
c multiply precondition matrix |
93 |
|
|
c------------------------------------ |
94 |
|
|
if (mupd .ne. 0) then |
95 |
|
|
call dostore( nn, xdiff, .false., 3 ) |
96 |
|
|
do i = 1, nn |
97 |
|
|
dd(i) = dd(i)*xdiff(i) |
98 |
|
|
end do |
99 |
|
|
end if |
100 |
|
|
|
101 |
|
|
c------------------------------------ |
102 |
|
|
c compute left hand side |
103 |
|
|
c------------------------------------ |
104 |
|
|
do j = jmin,jfin |
105 |
|
|
jp = j |
106 |
|
|
if (jp .gt. mupd) jp = jp-mupd |
107 |
|
|
call dostore( nn, xdiff, .false., 2*jp+2 ) |
108 |
|
|
|
109 |
|
|
temp = DDOT( nn, dd, 1, xdiff, 1 ) |
110 |
|
|
call MPI_ALLREDUCE(temp,r,1,MPI_DOUBLE_PRECISION,MPI_SUM, |
111 |
|
|
& MPI_COMM_WORLD,ierr) |
112 |
|
|
r = alpha(jp) - r |
113 |
|
|
! r = alpha(jp) - DDOT( nn, dd,1 , xdiff, 1 ) |
114 |
|
|
|
115 |
|
|
call dostore( nn, xdiff, .false., 2*jp+3 ) |
116 |
|
|
do i = 1, nn |
117 |
|
|
dd(i) = dd(i) + r*xdiff(i) |
118 |
|
|
end do |
119 |
|
|
end do |
120 |
|
|
|
121 |
|
|
return |
122 |
|
|
|
123 |
|
|
end |