10 SUBROUTINE bicub_derivative ( ya, x1a, x2a, y1a, y2a, y12a, i_max, k_max )
16 REAL,
INTENT(IN) :: ya(r_grid,ndft)
17 REAL,
INTENT(IN) :: x1a(r_grid)
18 REAL,
INTENT(IN) :: x2a(ndft)
19 REAL,
INTENT(OUT) :: y1a(r_grid,ndft)
20 REAL,
INTENT(OUT) :: y2a(r_grid,ndft)
21 REAL,
INTENT(OUT) :: y12a(r_grid,ndft)
22 INTEGER,
INTENT(IN) :: i_max
23 INTEGER,
INTENT(IN) :: k_max
32 y1a(i,k) = (ya(i+1,k)-ya(i-1,k)) / (x1a(i+1)-x1a(i-1))
33 y2a(i,k) = (ya(i,k+1)-ya(i,k-1)) / (x2a(k+1)-x2a(k-1))
34 y12a(i,k)= (ya(i+1,k+1)-ya(i+1,k-1)-ya(i-1,k+1)+ya(i-1,k-1)) &
35 /((x1a(i+1)-x1a(i-1))*(x2a(k+1)-x2a(k-1)))
56 y1a(i,k) = (ya(i,k)-ya(i-1,k)) / (x1a(i)-x1a(i-1))
57 y2a(i,k) = (ya(i,k+1)-ya(i,k-1)) / (x2a(k+1)-x2a(k-1))
58 y12a(i,k)= (ya(i,k+1)-ya(i,k-1)-ya(i-1,k+1)+ya(i-1,k-1)) &
59 /((x1a(i)-x1a(i-1))*(x2a(k+1)-x2a(k-1)))
65 y1a(i,k) = (ya(i+1,k)-ya(i-1,k)) / (x1a(i+1)-x1a(i-1))
66 y2a(i,k) = (ya(i,k)-ya(i,k-1)) / (x2a(k)-x2a(k-1))
67 y12a(i,k)= (ya(i+1,k)-ya(i+1,k-1)-ya(i-1,k)+ya(i-1,k-1)) &
68 /((x1a(i+1)-x1a(i-1))*(x2a(k)-x2a(k-1)))
73 y1a(i,k) = (ya(i,k)-ya(i-1,k)) / (x1a(i)-x1a(i-1))
74 y2a(i,k) = (ya(i,k)-ya(i,k-1)) / (x2a(k)-x2a(k-1))
75 y12a(i,k)= (ya(i,k)-ya(i,k-1)-ya(i-1,k)+ya(i-1,k-1)) &
76 /((x1a(i)-x1a(i-1))*(x2a(k)-x2a(k-1)))
78 END SUBROUTINE bicub_derivative
85 SUBROUTINE bicub_c ( ya, x1a, x2a, y1a, y2a, y12a, c_bicub, i_max, k_max )
91 REAL,
INTENT(IN) :: ya(r_grid,ndft)
92 REAL,
INTENT(IN) :: x1a(r_grid)
93 REAL,
INTENT(IN) :: x2a(ndft)
94 REAL,
INTENT(IN) :: y1a(r_grid,ndft)
95 REAL,
INTENT(IN) :: y2a(r_grid,ndft)
96 REAL,
INTENT(IN) :: y12a(r_grid,ndft)
97 REAL,
INTENT(OUT) :: c_bicub(r_grid,ndft,4,4)
98 INTEGER,
INTENT(IN) :: i_max
99 INTEGER,
INTENT(IN) :: k_max
102 INTEGER :: i, k, m, n
103 REAL :: y(4),y1(4),y2(4),y12(4),x1l,x1u,x2l,x2u
134 CALL bcucof(y,y1,y2,y12,x1u-x1l,x2u-x2l,c)
137 c_bicub(i,k,m,n)=c(m,n)
144 END SUBROUTINE bicub_c
150 SUBROUTINE bcucof(y,y1,y2,y12,d1,d2,c)
155 REAL,
INTENT(IN) :: y(4)
156 REAL,
INTENT(IN) :: y1(4)
157 REAL,
INTENT(IN) :: y2(4)
158 REAL,
INTENT(IN) :: y12(4)
159 REAL,
INTENT(IN) :: d1
160 REAL,
INTENT(IN) :: d2
161 REAL,
INTENT(OUT) :: c(4,4)
165 REAL :: d1d2,xx,cl(16),wt(16,16),x(16)
167 DATA wt/1,0,-3,2,4*0,-3,0,9,-6,2,0,-6,4,8*0,3,0,-9,6,-2,0,6,-4,10* &
168 0,9,-6,2*0,-6,4,2*0,3,-2,6*0,-9,6,2*0,6,-4,4*0,1,0,-3,2,-2,0,6,-4, &
169 1,0,-3,2,8*0,-1,0,3,-2,1,0,-3,2,10*0,-3,2,2*0,3,-2,6*0,3,-2,2*0, &
170 -6,4,2*0,3,-2,0,1,-2,1,5*0,-3,6,-3,0,2,-4,2,9*0,3,-6,3,0,-2,4,-2, &
171 10*0,-3,3,2*0,2,-2,2*0,-1,1,6*0,3,-3,2*0,-2,2,5*0,1,-2,1,0,-2,4, &
172 -2,0,1,-2,1,9*0,-1,2,-1,0,1,-2,1,10*0,1,-1,2*0,-1,1,6*0,-1,1,2*0, &
181 x(i+12) = y12(i)*d1d2
186 xx = xx + wt(i,k) * x(k)
198 END SUBROUTINE bcucof
205 SUBROUTINE bi_cub_spline ( rho_rdf, xg, ya, x1a, x2a, y1a, y2a, y12a, &
206 c_bicub, rdf, dg_drho, dg_dr, i_max, ih, k )
212 REAL,
INTENT(IN OUT) :: rho_rdf
213 REAL,
INTENT(IN OUT) :: xg
214 REAL,
INTENT(IN) :: ya(r_grid,ndft)
215 REAL,
INTENT(IN) :: x1a(r_grid)
216 REAL,
INTENT(IN) :: x2a(ndft)
217 REAL,
INTENT(IN) :: y1a(r_grid,ndft)
218 REAL,
INTENT(IN) :: y2a(r_grid,ndft)
219 REAL,
INTENT(IN) :: y12a(r_grid,ndft)
220 REAL,
INTENT(IN) :: c_bicub(r_grid,ndft,4,4)
221 REAL,
INTENT(OUT) :: rdf
222 REAL,
INTENT(OUT) :: dg_drho
223 REAL,
INTENT(OUT) :: dg_dr
224 INTEGER,
INTENT(IN OUT) :: i_max
226 INTEGER,
INTENT(OUT) :: ih
227 INTEGER,
INTENT(IN) :: k
232 REAL :: y(4),y1(4),y2(4),y12(4),x1l,x1u,x2l,x2u
236 IF ( rho_rdf < x1a(1) )
THEN 242 IF ( x1a(ih) <= rho_rdf .AND. rho_rdf < x1a(ih+1) )
GO TO 10
244 IF ( x1a(ih-1) <= rho_rdf .AND. rho_rdf < x1a(ih) )
THEN 250 CALL hunt ( x1a, i_max, rho_rdf, ih )
253 IF ( x2a(k) /= xg )
THEN 270 y1(3) = y1a(ih+1,k+1)
275 y2(3) = y2a(ih+1,k+1)
279 y12(2) = y12a(ih+1,k)
280 y12(3) = y12a(ih+1,k+1)
281 y12(4) = y12a(ih,k+1)
290 c(m,n) = c_bicub( ih, k, m, n )
293 CALL bcuint ( x1l, x1u, x2l, x2u, rho_rdf, xg, c, rdf, dg_drho, dg_dr )
295 END SUBROUTINE bi_cub_spline
307 SUBROUTINE hunt ( xx, n, x, jlo )
312 INTEGER,
INTENT(IN) :: n
313 INTEGER,
INTENT(OUT) :: jlo
314 REAL,
INTENT(IN) :: xx(n)
318 INTEGER :: inc,jhi,jm
322 ascnd = xx(n) >= xx(1)
323 IF( jlo <= 0 .OR. jlo > n )
THEN 329 IF( x >= xx(jlo) .EQV. ascnd )
THEN 333 ELSE IF ( x >= xx(jhi) .EQV. ascnd )
THEN 343 ELSE IF ( x < xx(jlo) .EQV. ascnd )
THEN 349 3
IF (jhi-jlo == 1 )
THEN 350 IF ( x == xx(n)) jlo = n - 1
351 IF ( x == xx(1) ) jlo = 1
354 jm = ( jhi + jlo ) / 2
355 IF ( x >= xx(jm) .EQV. ascnd )
THEN 370 SUBROUTINE bcuint ( x1l, x1u, x2l, x2u, x1, x2, c, ansy, ansy1, ansy2 )
379 REAL,
INTENT(IN OUT) :: x1l
380 REAL,
INTENT(IN OUT) :: x1u
381 REAL,
INTENT(IN OUT) :: x2l
382 REAL,
INTENT(IN OUT) :: x2u
383 REAL,
INTENT(IN OUT) :: x1
384 REAL,
INTENT(IN OUT) :: x2
385 REAL,
INTENT(IN) :: c(4,4)
386 REAL,
INTENT(OUT) :: ansy
387 REAL,
INTENT(OUT) :: ansy1
388 REAL,
INTENT(OUT) :: ansy2
398 IF ( x1u == x1l .OR. x2u == x2l ) pause
'bad input in bcuint' 399 t = (x1-x1l) / (x1u-x1l)
400 u = (x2-x2l) / (x2u-x2l)
405 ansy = t *ansy + ( (c(i,4)*u + c(i,3))*u+c(i,2) )*u + c(i,1)
406 ansy2 = t *ansy2 + ( 3.*c(i,4)*u+2.*c(i,3) )*u + c(i,2)
407 ansy1 = u *ansy1 + ( 3.*c(4,i)*t+2.*c(3,i) )*t + c(2,i)
409 ansy1 = ansy1 / (x1u-x1l)
410 ansy2 = ansy2 / (x2u-x2l)
412 END SUBROUTINE bcuint
430 SUBROUTINE spline ( x, y, n, yp1, ypn, y2 )
435 INTEGER,
INTENT(IN) :: n
436 REAL,
INTENT(IN) :: x(n)
437 REAL,
INTENT(IN) :: y(n)
438 REAL,
INTENT(IN) :: yp1
439 REAL,
INTENT(IN) :: ypn
440 REAL,
INTENT(OUT) :: y2(n)
443 INTEGER,
PARAMETER :: nmax = 1000
445 REAL :: p, qn,
sig, un, u(nmax)
448 If(n > nmax) stop
'Increase NMAX in Spline!!' 452 IF ( yp1 > 0.99e30 )
THEN 457 u(1) = ( 3.0/(x(2)-x(1)) ) * ( (y(2)-y(1))/(x(2)-x(1))-yp1 )
460 IF ( (x(i+1)-x(i)) == 0.0 .OR. (x(i)-x(i-1)) == 0.0 .OR. (x(i+1)-x(i-1)) == 0.0 )
THEN 461 write (*,*)
'Surface Tension Code: error in spline-interpolation' 464 sig = (x(i)-x(i-1)) / (x(i+1)-x(i-1))
465 p =
sig*y2(i-1) + 2.0
466 y2(i) = (
sig-1.0) / p
467 u(i) = ( 6.0 * ((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1))) / (x(i+1)-x(i-1)) &
470 IF ( ypn > 0.99e30 )
THEN 475 un = (3.0/(x(n)-x(n-1))) * (ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
477 y2(n) = (un-qn*u(n-1)) / (qn*y2(n-1)+1.0)
479 y2(k) = y2(k) * y2(k+1) + u(k)
482 END SUBROUTINE spline
493 SUBROUTINE splint_integral ( xa, ya, y2a, n, xlo, xhi, integral )
498 INTEGER,
INTENT(IN) :: n
499 REAL,
INTENT(IN) :: xa(n)
500 REAL,
INTENT(IN) :: ya(n)
501 REAL,
INTENT(IN) :: y2a(n)
502 REAL,
INTENT(IN) :: xlo
503 REAL,
INTENT(IN) :: xhi
504 REAL,
INTENT(OUT) :: integral
507 INTEGER :: k, khi_l, klo_l, khi_h, klo_h
508 REAL :: xl, xh, h, int, x0, x1, y0, y1, y20, y21
514 1
IF ( khi_l-klo_l > 1 )
THEN 515 k = ( khi_l + klo_l ) / 2
516 IF ( xa(k) > xlo )
THEN 526 2
IF ( khi_h-klo_h > 1 )
THEN 527 k = ( khi_h + klo_h ) / 2
528 IF ( xa(k) > xhi )
THEN 542 IF ( khi_h > khi_l )
THEN 544 ELSE IF ( khi_h == khi_l )
THEN 547 WRITE (*,*)
'error in spline-integration' 551 h = xa(khi_l) - xa(klo_l)
552 IF ( h == 0.0 ) pause
'bad xa input in splint' 569 int = -1.0/h * ( xl*((x1-.5*xl)*y0 + (0.5*xl-x0)*y1) &
570 -y20/24.*(x1-xl)**4 + y20/6.*(0.5*xl*xl-x1*xl)*h*h &
571 +y21/24.*(xl-x0)**4 - y21/6.*(0.5*xl*xl-x0*xl)*h*h )
572 int = int + 1.0/h * ( xh*((x1-.5*xh)*y0 + (0.5*xh-x0)*y1) &
573 -y20/24.*(x1-xh)**4 + y20/6.*(0.5*xh*xh-x1*xh)*h*h &
574 +y21/24.*(xh-x0)**4 - y21/6.*(0.5*xh*xh-x0*xh)*h*h )
576 integral = integral + int
581 IF (khi_h /= (khi_l-1))
GO TO 3
583 END SUBROUTINE splint_integral
592 REAL FUNCTION praxis( t0, machep, h0, n, prin, x, f, fmin )
597 REAL,
INTENT(IN OUT) :: t0
598 REAL,
INTENT(IN) :: machep
599 REAL,
INTENT(IN) :: h0
601 INTEGER,
INTENT(IN OUT) :: prin
602 REAL,
INTENT(IN OUT) :: x(n)
604 REAL,
INTENT(IN OUT) :: fmin
661 INTEGER :: nl,nf,kl,kt,ktm,idim,i,j,k,k2,km1,klmk,ii,im1
662 REAL :: s,sl,dn,dmin,fx,f1,lds,ldt,t,h,sf,df,qf1,qd0, qd1,qa,qb,qc
663 REAL :: m2,m4,small,vsmall,large,vlarge,scbd,ldfac,t2, dni,value
670 REAL :: d(20),y(20),z(20),q0(20),q1(20),v(20,20)
671 COMMON /global/ fx,ldt,dmin,nf,nl /q/ v,q0,q1,qa,qb,qc,qd0,qd1,qf1
683 small = machep*machep
704 IF (illc) ldfac = 0.1
714 IF (h < 100*t) h = 100*t
729 IF (prin > 0)
CALL print(n,x,prin,fmin)
739 CALL min(n,1,2,d(1),s,
value,.false.,f,x,t,machep,h)
740 IF (s > 0.d0)
GO TO 50
744 50
IF (sf > 0.9d0*d(1).AND.0.9d0*sf < d(1))
GO TO 70
755 IF (kt > 0) illc=.true.
763 IF(.NOT.illc)
GO TO 95
765 s=(0.1d0*ldt+t2*(10**kt))*(random(n)-0.5d0)
780 CALL min(n,k2,2,d(k2),s,
value,.false.,f,x,t,machep,h)
784 97 s=d(k2)*((s+z(k2))**2)
789 IF (illc.OR.(df >= abs((100*machep)*fx)))
GO TO 110
796 110
IF (k == 2.AND.prin > 1)
CALL vcprnt(1,d,n)
804 CALL min(n,k2,2,d(k2),s,
value,.false.,f,x,t,machep,h)
817 IF (lds <= small)
GO TO 160
824 IF (klmk < 1)
GO TO 141
841 CALL min(n,k,4,d(k),lds,
value,.true.,f,x,t,machep,h)
842 IF (lds > 0.d0)
GO TO 160
848 IF (ldt < lds) ldt=lds
849 IF (prin > 0)
CALL print(n,x,prin,fmin)
859 IF (ldt > (0.5*t2)) kt=-1
861 IF (kt > ktm)
GO TO 400
867 CALL quad(n,f,x,t,machep,h)
871 IF (dn < d(i)) dn=d(i)
873 IF (prin > 3)
CALL maprnt(1,v,idim,n)
883 IF (scbd <= 1.d0)
GO TO 200
891 IF (z(i) < m4) z(i)=m4
897 IF (z(i) <= scbd)
GO TO 189
923 CALL minfit(idim,n,machep,vsmall,v,d)
927 IF (scbd <= 1.d0)
GO TO 250
949 IF (dni > large)
GO TO 265
950 IF (dni < small)
GO TO 260
960 CALL sort(idim,n,d,v)
962 IF (dmin < small) dmin=small
964 IF (m2*d(1) > dmin) illc=.true.
965 IF (prin > 1.AND.scbd > 1.d0)
CALL vcprnt(2,z,n)
966 IF (prin > 1)
CALL vcprnt(3,d,n)
967 IF (prin > 3)
CALL maprnt(2,v,idim,n)
974 400
IF (prin > 0)
CALL vcprnt(4,x,n)
984 SUBROUTINE minfit(m,n,machep,tol,ab,q)
988 INTEGER,
INTENT(IN OUT) :: m
989 INTEGER,
INTENT(IN) :: n
990 REAL,
INTENT(IN) :: machep
991 REAL,
INTENT(IN OUT) :: tol
992 REAL,
INTENT(IN OUT) :: ab(m,n)
993 REAL,
INTENT(OUT) :: q(n)
994 INTEGER :: i,j,k,l, kk,kt,l2,ll2,ii,lp1
998 REAL :: x,eps,e(20),g,s, f,h,y,c,z,temp
1006 IF (n == 1)
GO TO 200
1018 IF (s < tol)
GO TO 4
1021 IF (f >= 0.d0) g = -g
1028 f = f + ab(k,i)*ab(k,j)
1032 ab(k,j) = ab(k,j) + f*ab(k,i)
1039 s = s + ab(i,j)*ab(i,j)
1042 IF (s < tol)
GO TO 10
1043 IF (i == n)
GO TO 16
1046 IF (f >= 0.d0) g = -g
1048 IF (i == n)
GO TO 10
1056 s = s + ab(j,k)*ab(i,k)
1059 ab(j,k) = ab(j,k) + s*e(k)
1062 10 y = abs(q(i)) + abs(e(i))
1071 IF (g == 0.d0)
GO TO 23
1079 s = s + ab(i,k)*ab(k,j)
1082 ab(k,j) = ab(k,j) + s*ab(k,i)
1099 IF (kt <= 30)
GO TO 102
1102 1000
FORMAT (
' QR FAILED')
1106 IF (abs(e(l)) <= eps)
GO TO 120
1108 IF (abs(q(l-1)) <= eps)
EXIT 1116 IF (abs(f) <= eps)
GO TO 120
1119 IF (abs(f) < abs(g))
GO TO 113
1127 112 h = abs(f)*sqrt(1 + (g/f)**2)
1129 113 h = abs(g)*sqrt(1 + (f/g)**2)
1131 IF (h /= 0.d0)
GO TO 115
1139 IF (l == k)
GO TO 140
1145 f = ((y - z)*(y + z) + (g - h)*(g + h))/(2*h*y)
1146 g = sqrt(f*f + 1.0d0)
1148 IF (f >= 0.d0) temp = f + g
1149 f = ((x - z)*(x + z) + h*(y/temp - h))/x
1154 IF (lp1 > k)
GO TO 133
1160 IF (abs(f) < abs(h))
GO TO 123
1168 122 z = abs(f)*sqrt(1 + (h/f)**2)
1170 123 z = abs(h)*sqrt(1 + (f/h)**2)
1172 IF (z /= 0.d0)
GO TO 125
1184 ab(j,i-1) = x*c + z*s
1185 ab(j,i) = -x*s + z*c
1187 IF (abs(f) < abs(h))
GO TO 129
1195 128 z = abs(f)*sqrt(1 + (h/f)**2)
1197 129 z = abs(h)*sqrt(1 + (f/h)**2)
1199 IF (z /= 0.d0)
GO TO 131
1212 140
IF (z >= 0.d0) cycle
1222 END SUBROUTINE minfit
1231 SUBROUTINE min(n,j,nits,d2,x1,f1,fk,f,x,t,machep,h)
1235 INTEGER,
INTENT(IN) :: n
1238 REAL,
INTENT(IN OUT) :: d2
1239 REAL,
INTENT(IN OUT) :: x1
1240 REAL,
INTENT(IN OUT) :: f1
1243 REAL,
INTENT(IN OUT) :: x(n)
1244 REAL,
INTENT(IN) :: t
1245 REAL,
INTENT(IN) :: machep
1246 REAL,
INTENT(IN) :: h
1253 REAL :: small,sf1,sx1,s,temp, xm,x2,f2,d1
1258 REAL :: v(20,20),q0(20),q1(20),qa,qb,qc,qd0,qd1,qf1
1259 COMMON /global/ fx,ldt,dmin,nf,nl /q/ v,q0,q1,qa,qb,qc,qd0,qd1,qf1
1293 t2 = m4*sqrt(abs(fx)/temp + s*ldt) + m2*ldt
1295 IF (dz.AND.t2 > s) t2 = s
1296 t2 = dmax1(t2,small)
1297 t2 = dmin1(t2,.01d0*h)
1298 IF (.NOT.fk.OR.f1 > fm)
GO TO 2
1301 2
IF (fk.AND.abs(x1) >= t2)
GO TO 3
1303 IF (x1 < 0.d0) temp=-1.d0
1305 f1 = flin(n,j,x1,f,x,nf)
1306 3
IF (f1 > fm)
GO TO 4
1309 4
IF (.NOT.dz)
GO TO 6
1312 IF (f0 >= f1) x2 = 2.d0*x1
1313 f2 = flin(n,j,x2,f,x,nf)
1314 IF (f2 > fm)
GO TO 5
1317 5 d2 = (x2*(f1 - f0)-x1*(f2 - f0))/((x1*x2)*(x1 - x2))
1319 6 d1 = (f1 - f0)/x1 - x1*d2
1322 IF (d2 > small)
GO TO 7
1324 IF (d1 >= 0.d0) x2 = -x2
1326 7 x2 = (-.5d0*d1)/d2
1327 8
IF (abs(x2) <= h)
GO TO 11
1335 11 f2 = flin(n,j,x2,f,x,nf)
1336 IF (k >= nits.OR.f2 <= f0)
GO TO 12
1339 IF (f0 < f1.AND.(x1*x2) > 0.d0)
GO TO 4
1344 IF (f2 <= fm)
GO TO 13
1349 14
IF (abs(x2*(x2 - x1)) <= small)
GO TO 15
1350 d2 = (x2*(f1-f0) - x1*(fm-f0))/((x1*x2)*(x1 - x2))
1352 15
IF (k > 0) d2 = 0.d0
1353 16
IF (d2 <= small) d2 = small
1356 IF (sf1 >= fx)
GO TO 17
1360 17
IF (j == 0)
RETURN 1362 x(i) = x(i) + x1*v(i,j)
1372 SUBROUTINE vcprnt(option,v,n)
1376 REAL,
INTENT(IN OUT) :: v(n)
1381 SELECT CASE ( option )
1392 1
WRITE (6,101) (v(i),i=1,n)
1394 2
WRITE (6,102) (v(i),i=1,n)
1396 3
WRITE (6,103) (v(i),i=1,n)
1398 4
WRITE (6,104) (v(i),i=1,n)
1400 101
FORMAT (/
' THE SECOND DIFFERENCE ARRAY D(*) IS:'/ (e32.14,4e25.14))
1401 102
FORMAT (/
' THE SCALE FACTORS ARE:'/(e32.14,4e25.14))
1402 103
FORMAT (/
' THE APPROXIMATING QUADR. FORM HAS PRINCIPAL VALUES:'/ &
1404 104
FORMAT (/
' X IS:',e26.14/(e32.14))
1405 END SUBROUTINE vcprnt
1412 SUBROUTINE print(n,x,prin,fmin)
1416 REAL,
INTENT(IN OUT) :: x(n)
1417 INTEGER,
INTENT(IN OUT) :: prin
1418 REAL,
INTENT(IN OUT) :: fmin
1425 COMMON /global/ fx,ldt,dmin,nf,nl
1427 WRITE (6,101) nl,nf,fx
1429 IF (fx <= fmin)
GO TO 1
1431 WRITE (6,102) fmin,ln
1433 1
WRITE (6,103) fmin
1434 2
IF (n > 4.AND.prin <= 2)
RETURN 1435 WRITE (6,104) (x(i),i=1,n)
1437 101
FORMAT (/
' AFTER',i6, &
1438 ' LINEAR SEARCHES, THE FUNCTION HAS BEEN EVALUATED',i6, &
1439 ' TIMES. THE SMALLEST VALUE FOUND IS F(X) = ',e21.14)
1440 102
FORMAT (
' LOG (F(X)-',e21.14,
') = ',e21.14)
1441 103
FORMAT (
' LOG (F(X)-',e21.14,
') IS UNDEFINED.')
1442 104
FORMAT (
' X IS:',e26.14/(e32.14))
1443 END SUBROUTINE print
1450 SUBROUTINE maprnt(option,v,m,n)
1454 REAL,
INTENT(IN OUT) :: v(m,n)
1455 INTEGER,
INTENT(IN OUT) :: m
1456 INTEGER,
INTENT(IN) :: n
1465 SELECT CASE ( option )
1472 101
FORMAT (/
' THE NEW DIRECTIONS ARE:')
1475 102
FORMAT (
' AND THE PRINCIPAL AXES:')
1476 3
IF (n < upp) upp = n
1478 WRITE (6,104) (v(i,j),j=low,upp)
1486 104
FORMAT (e32.14,4e25.14)
1487 END SUBROUTINE maprnt
1494 REAL FUNCTION random(naught)
1497 INTEGER,
INTENT(IN OUT) :: naught
1499 REAL :: ran1,ran3(127),half
1500 INTEGER :: i,j,ran2,q,r
1503 SAVE init,ran2,ran1,ran3
1506 r = mod(naught,8190) + 1
1512 r = mod(1756*r,8191)
1514 ran1 = (ran1 + q)*(1.0d0/256)
1519 3
IF (ran2 == 1) ran2 = 128
1521 ran1 = ran1 + ran3(ran2)
1523 IF (ran1 >= 0.d0) half = -half
1526 random = ran1 + .5d0
1536 REAL FUNCTION flin (n,j,l,f,x,nf)
1539 INTEGER,
INTENT(IN) :: n
1540 INTEGER,
INTENT(IN OUT) :: j
1541 REAL,
INTENT(IN) :: l
1543 REAL,
INTENT(IN) :: x(n)
1544 INTEGER,
INTENT(OUT) :: nf
1554 REAL :: v(20,20),q0(20),q1(20),qa,qb,qc,qd0,qd1,qf1
1555 COMMON /q/ v,q0,q1,qa,qb,qc,qd0,qd1,qf1
1560 t(i) = x(i) + l*v(i,j)
1564 2 qa = (l*(l - qd1))/(qd0*(qd0 + qd1))
1565 qb = ((l + qd0)*(qd1 - l))/(qd0*qd1)
1566 qc = (l*(l + qd0))/(qd1*(qd0 + qd1))
1568 t(i) = (qa*q0(i) + qb*x(i)) + qc*q1(i)
1581 SUBROUTINE sort(m,n,d,v)
1584 INTEGER,
INTENT(IN OUT) :: m
1585 INTEGER,
INTENT(IN) :: n
1586 REAL,
INTENT(IN OUT) :: d(n)
1587 REAL,
INTENT(IN OUT) :: v(m,n)
1589 INTEGER :: i,j,k,nm1,ip1
1601 IF (d(j) <= s) cycle
1621 SUBROUTINE quad(n,f,x,t,machep,h)
1624 INTEGER,
INTENT(IN) :: n
1626 REAL,
INTENT(IN OUT) :: x(n)
1627 REAL,
INTENT(IN OUT) :: t
1629 REAL,
INTENT(IN OUT) :: h
1641 REAL :: v(20,20),q0(20),q1(20),qa,qb,qc,qd0,qd1,qf1
1642 COMMON /global/ fx,ldt,dmin,nf,nl /q/ v,q0,q1,qa,qb,qc,qd0,qd1,qf1
1653 qd1 = qd1 + (s-l)**2
1658 IF (qd0 <= 0.d0 .OR. qd1 <= 0.d0 .OR. nl < 3*n*n)
GO TO 2
1660 CALL min(n,0,2,s,l,
value,.true.,f,x,t,machep,h)
1661 qa = (l*(l-qd1))/(qd0*(qd0+qd1))
1662 qb = ((l+qd0)*(qd1-l))/(qd0*qd1)
1663 qc = (l*(l+qd0))/(qd1*(qd0+qd1))
1673 x(i) = (qa*s + qb*x(i)) + qc*q1(i)
double sig
correlated to variance of initial distribution
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW Module DFT_MODULE This module...