1*
2*     $Id$
3*
4
5      subroutine paw_lmbfgs_init(max_m0,g)
6      implicit none
7      integer max_m0
8      complex*16 g(*)
9
10#include "paw_lmbfgs_common.fh"
11
12
13*     **** external functions ****
14      integer  paw_psi_ne
15      external paw_psi_ne
16
17      max_m = max_m0
18      call Pack_npack(1,npack1)
19      neall = paw_psi_ne(1)+paw_psi_ne(2)
20      nsize = neall*npack1
21      call paw_Grsm_list_init('lmbfgs',2*max_m,2*nsize)
22
23      m = 1
24c      call Grsm_gg_dScale(npack1,neall,(-1.0d0),g,g)
25      call Grsm_gg_dScale1(npack1,neall,(-1.0d0),g)
26      call paw_Grsm_list_store(2*m-1,g)
27      call paw_Grsm_list_store(2*m,g)
28c      call Grsm_gg_dScale(npack1,neall,(-1.0d0),g,g)
29      call Grsm_gg_dScale1(npack1,neall,(-1.0d0),g)
30      return
31      end
32
33
34
35      subroutine paw_lmbfgs(tmin,g,s)
36      implicit none
37      real*8 tmin
38      complex*16 g(*)
39      complex*16 s(*)
40
41
42#include "bafdecls.fh"
43#include "paw_lmbfgs_common.fh"
44
45*     **** local variables ****
46      logical value,ondisk
47      integer k,ngras
48      real*8 alpha(20),beta,sum,h0
49      integer yy(2),ss(2)
50
51*     **** external functions ****
52      logical  control_lmbfgs_ondisk
53      external control_lmbfgs_ondisk
54
55      ondisk = control_lmbfgs_ondisk()
56      ngras = 2*npack1*neall
57
58      !call Grsm_gg_dScale(npack1,neall,(-1.0d0),g,g)
59      call dscal(ngras,(-1.0d0),g,1)
60
61*     **** allocate yy and ss ****
62      if (ondisk) then
63        value = BA_push_get(mt_dbl,2*nsize,'yy',yy(2),yy(1))
64        value = value.and.
65     >          BA_push_get(mt_dbl,2*nsize,'ss',ss(2),ss(1))
66        if (.not.value) call errquit('paw_lmbfgs: push stack',0,0)
67      end if
68      call dcopy(ngras,g,1,s,1)
69
70      if (ondisk) then
71        call paw_Grsm_list_load(2*m-1,dbl_mb(yy(1)))
72        call paw_Grsm_list_load(2*m,  dbl_mb(ss(1)))
73      else
74        call paw_Grsm_list_ptr(2*m-1,yy(1))
75        call paw_Grsm_list_ptr(2*m,  ss(1))
76      end if
77      call paw_psi_1geodesic_Gtransport(tmin,dbl_mb(yy(1)))
78      call paw_psi_1geodesic_transport(tmin,dbl_mb(ss(1)))
79      call daxpy(ngras,(-1.0d0),g,1,dbl_mb(yy(1)),1)
80      call dscal(ngras,(-1.0d0),dbl_mb(yy(1)),1)
81      if (ondisk) then
82        call paw_Grsm_list_store(2*m-1,dbl_mb(yy(1)))
83        call paw_Grsm_list_store(2*m,  dbl_mb(ss(1)))
84      end if
85
86      call Grsm_gg_trace(npack1,neall,
87     >                   dbl_mb(yy(1)),
88     >                   dbl_mb(ss(1)),
89     >                   sum)
90      !*** exit if dividing by small number ***
91      if (dabs(sum).lt.1.0d-15) go to 211
92
93      rho(m) = 1.0d0/sum
94      call Grsm_gg_trace(npack1,neall,
95     >                   dbl_mb(ss(1)),
96     >                   s,
97     >                   sum)
98      alpha(m) = rho(m)*sum
99      call daxpy(ngras,(-alpha(m)),dbl_mb(yy(1)),1,s,1)
100      do k=(m-1),1, -1
101         if (ondisk) then
102           call paw_Grsm_list_load(2*k-1,dbl_mb(yy(1)))
103           call paw_Grsm_list_load(2*k,  dbl_mb(ss(1)))
104         else
105           call paw_Grsm_list_ptr(2*k-1,yy(1))
106           call paw_Grsm_list_ptr(2*k,  ss(1))
107         end if
108         call paw_psi_1geodesic_Gtransport(tmin,dbl_mb(yy(1)))
109         call paw_psi_1geodesic_Gtransport(tmin,dbl_mb(ss(1)))
110         if (ondisk) then
111           call paw_Grsm_list_store(2*k-1,dbl_mb(yy(1)))
112           call paw_Grsm_list_store(2*k,  dbl_mb(ss(1)))
113         end if
114
115         call Grsm_gg_trace(npack1,neall,
116     >                   dbl_mb(ss(1)),
117     >                   s,
118     >                   sum)
119         alpha(k) = rho(k)*sum
120         call daxpy(ngras,(-alpha(k)),dbl_mb(yy(1)),1,s,1)
121      end do
122
123
124*     **** preconditioner ****
125c      call Grsm_gg_dScale(npack1,neall,h0,s,s)
126
127
128      do k=1,(m-1)
129
130         call Grsm_gg_trace(npack1,neall,
131     >                   dbl_mb(yy(1)),
132     >                   s,
133     >                   sum)
134         beta = rho(k)*sum
135         sum = alpha(k) - beta
136         call daxpy(ngras,sum,dbl_mb(ss(1)),1,s,1)
137
138         if (ondisk) then
139           call paw_Grsm_list_load(2*(k+1)-1,dbl_mb(yy(1)))
140           call paw_Grsm_list_load(2*(k+1),  dbl_mb(ss(1)))
141         else
142           call paw_Grsm_list_ptr(2*(k+1)-1,yy(1))
143           call paw_Grsm_list_ptr(2*(k+1),  ss(1))
144         end if
145
146      end do
147      call Grsm_gg_trace(npack1,neall,
148     >                dbl_mb(yy(1)),
149     >                s,
150     >                sum)
151      beta = rho(m)*sum
152      sum = alpha(m) - beta
153      call daxpy(ngras,sum,dbl_mb(ss(1)),1,s,1)
154      if (m.lt.max_m) then
155         m = m+1
156      else
157         call paw_Grsm_list_shift()
158         call paw_Grsm_list_shift()
159         do k=1,(m-1)
160            rho(k) = rho(k+1)
161         end do
162      end if
163
164      call dscal(ngras,(-1.0d0),s,1)
165 211  call paw_Grsm_list_store(2*m-1,g)
166      call paw_Grsm_list_store(2*m,s)
167      call dscal(ngras, (-1.0d0),g,1)
168
169      if (ondisk) then
170        value = BA_pop_stack(ss(2))
171        value = value.and.
172     >          BA_pop_stack(yy(2))
173        if (.not.value) call errquit('paw_lmbfgs:pop stack',2,0)
174      end if
175
176      return
177      end
178
179