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