1 double precision function linesrch(n,x,f,g) 2* 3* $Id$ 4* 5 implicit none 6 integer n 7 double precision x(n) 8 double precision f(n) 9 double precision g(n) 10 integer ia, ib, ic 11 double precision dx, tol, prec 12 double precision step 13 logical linebracket 14 external linebracket 15 16 write(6,*) 17 write(6,*) 18 do ia=max(1,(n-4)),n 19 write(6,899) ia,x(ia),f(ia),g(ia) 20 899 format(i5,2f24.12,e20.6) 21 enddo 22 if (n.eq.2) then 23 dx = (g(n) - g(n-1))/(x(n)-x(n-1)) 24 write(6,884) dx 25 884 format('Curvature:',e15.4) 26 if ((dx.lt.0.d0).and.(f(n).gt.f(n-1))) then 27 linesrch = -x(n) 28 return 29 endif 30 endif 31 tol = 1.d-13 32 step = (g(n)*x(n-1) - g(n-1)*x(n))/(g(n) - g(n-1)) 33c 34c Bracketing 35c 36 if (n.gt.2) then 37 if (linebracket(n,x,f,ia,ib,ic)) then 38 write(6,900) ib,x(ib),ia,x(ia),(f(ia)-f(ib)), 39 $ ic,x(ic),(f(ic)-f(ib)) 40 900 format( 'Minimum:',10x,i5,f22.14 41 $ /,'Left Bracket:',5x,i5,f22.14,5x,e14.6, 42 $ /,'Right Bracket:',4x,i5,f22.14,5x,e14.6) 43 write(6,771) step,n,n-1 44 771 format(5x,'Step:',f20.10,' from ',i3,'-',i3) 45c 46c If last step is not the minimum cannot use it 47c in NR extrapolation, use current brackets and 48c sign of gradient for extrapolation. 49c 50 if (ib.ne.n) then 51 if (g(ib).gt.0) then 52 step = (g(ia)*x(ib) - g(ib)*x(ia))/(g(ia) - g(ib)) 53 write(6,771) step,ia,ib 54 else 55 step = (g(ic)*x(ib) - g(ib)*x(ic))/(g(ic) - g(ib)) 56 write(6,771) step,ib,ic 57 endif 58 endif 59c 60c Apply bracketing criterion only if brackets 61c are greater than sqrt(relative machine-precision)*x 62c and difference in function values is greater than machine precision 63c 64 prec = abs(sqrt(tol)*x(ib)) 65 if (((step.lt.x(ia)).and.(abs(x(ib)-x(ia)).gt.prec).and. 66 $ (abs(f(ia)-f(ib)).gt.(tol*f(ib)))).or. 67 $ ((step.gt.x(ic)).and.(abs(x(ic)-x(ia)).gt.prec).and. 68 $ (abs(f(ic)-f(ib)).gt.(tol*f(ib))))) then 69 if (g(ib).gt.0) then 70 step = (x(ia)+x(ib))/2.d0 71 else 72 step = (x(ib)+x(ic))/2.d0 73 endif 74 write(6,772) step 75 772 format('Step outside bracket -- bisection step',f20.10) 76 endif 77 else 78c 79c Cannot bracket, assume function is monotonic 80c and use points with least function value 81c 82 write(6,774) ib,x(ib),f(ib) 83 774 format('Minimum: ',i5,2f20.10) 84 if ((ia.ne.0).and.(g(ib).lt.0.0d0)) then 85 dx = (g(ia) - g(ib))/(x(ia)-x(ib)) 86 write(6,775) ia,x(ia),f(ia),dx 87 775 format('Extrapolate right:',i5,2f20.10,5x, 88 $ '2nd deriv',2x,e12.3) 89 step = x(ib) + (x(ib) - x(1)) 90 else if ((ic.ne.0).and.(g(ib).gt.0.d0)) then 91 dx = (g(ic) - g(ib))/(x(ic)-x(ib)) 92 write(6,776) ic,x(ic),f(ic),dx 93 776 format('Extrapolate left:',i5,2f20.10,5x, 94 $ '2nd deriv',2x,e12.3) 95 step = x(ib) - (x(1) - x(ib)) 96 else 97 write(6,773) 98 773 format('Internal error....') 99 endif 100 endif 101 else if (n.eq.1) then 102 step = 1.d0 103 endif 104 linesrch = step 105 write(6,881) 106 881 format(//) 107 return 108 end 109 110 111 112 113 114 115 logical function linebracket(n,x,f,ia,ib,ic) 116 implicit none 117 integer n 118 double precision x(n),f(n) 119 integer ia,ib,ic 120 double precision xmin, ymin, del 121 integer i 122 123 ia = 0 124 ic = 0 125 ib = 1 126 xmin = f(1) 127 do i=2,n 128 if (f(i).lt.xmin) then 129 xmin = f(i) 130 ib = i 131 endif 132 enddo 133 xmin = -1.d0*abs(xmin) 134 ymin = abs(xmin) 135 do i=1,n 136 del = x(i) - x(ib) 137 if ((del.gt.0.d0).and.(del.lt.ymin)) then 138 ymin = del 139 ic = i 140 else if ((del.lt.0.d0).and.(del.gt.xmin)) then 141 xmin = del 142 ia = i 143 endif 144 enddo 145 linebracket = ((ia.ne.0).and.(ic.ne.0)) 146 return 147 end 148 149 150 151 152 153 154 155 156 157 158 double precision function linesrch1(n,x,f,g) 159 implicit none 160#include "global.fh" 161 integer n 162 double precision x(n) 163 double precision f(n) 164 double precision g(n) 165 double precision a, b1, b2, c1, c2, f1, f2, x0 166 double precision e1, e2, err 167 save e1, e2 168 data err/0.d0/ 169 170 if (n.eq.1) then 171 linesrch1 = 1.d0 172 return 173 else 174 if (n.gt.2) then 175 err = max(abs((f(n)-e1)/f(n)),abs((f(n)-e2)/f(n))) 176C if (ga_nodeid().eq.0) write(6,901) x(n),e1,e2,f(n),err 177 901 format(//,'Predictions @ ',f12.6,5x,2f20.10, 178 $ /,'Obs.',f20.10, 179 $ /,'Relative error:',e12.2) 180 endif 181 a = (g(n)-g(n-1))/(2.d0*(x(n)-x(n-1))) 182 b1 = g(n-1) - 4.d0*a*x(n-1) 183 b2 = g(n) - 4.d0*a*x(n) 184 c1 = f(n-1) - a*x(n-1)**2 - b1*x(n-1) 185 c2 = f(n) - a*x(n)**2 - b2*x(n) 186 f1 = a*x(n-1)**2 + b1*x(n-1) + c1 187 f2 = a*x(n)**2 + b2*x(n) + c2 188C if (ga_nodeid().eq.0) 189C $ write(6,900) a,b1,b2,c1,c2,f1,f(n-1),f2,f(n) 190 900 format('A: ',5x,f20.10,/, 191 $ 'B(1):',5x,f20.10,/, 192 $ 'B(2):',5x,f20.10,/, 193 $ 'C(1):',5x,f20.10,/, 194 $ 'C(2):',5x,f20.10,/, 195 $ 'F(1):',5x,f20.10,' f(1):',f20.10,/, 196 $ 'F(2):',5x,f20.10,' f(2):',f20.10) 197C if (err.lt.1.d-1) then 198 x0 = (g(n)*x(n-1) - g(n-1)*x(n))/(g(n) - g(n-1)) 199 e1 = a*x0**2 + b1*x0 + c1 200 e2 = a*x0**2 + b2*x0 + c2 201 linesrch1 = x0 202C else 203 204 endif 205 return 206 end 207 208 209 210