subroutine aptripq (px, py, pz,
& ac, ax, ay, az, axy, ayz, azx, axx, ayy, azz,
& bc, bx, by, bz, bxy, byz, bzx, bxx, byy, bzz,
& cc, cx, cy, cz, cxy, cyz, czx, cxx, cyy, czz,
& tol, tx, ty, tz, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTRIPQ
c
c call aptripq (px, py, pz,
c & ac, ax, ay, az, axy, ayz, azx, axx, ayy, azz,
c & bc, bx, by, bz, bxy, byz, bzx, bxx, byy, bzz,
c & cc, cx, cy, cz, cxy, cyz, czx, cxx, cyy, czz,
c & tol, tx, ty, tz, nerr)
c
c Version: aptripq Updated 1999 August 11 14:00.
c aptripq Originated 1997 December 10 15:40.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find, near point p = (px, py, pz), one of the 8 possible
c triple points of intersection t = (tx, ty, tz) of the three
c quadric surfaces with the implicit equations:
c
c A = ac + ax * x + ay * y + az * z
c + axy * x * y + ayz * y * z + azx * z * x
c + axx * x**2 + ayy * y**2 + azz * z**2 = 0
c B = bc + bx * x + by * y + bz * z
c + bxy * x * y + byz * y * z + bzx * z * x
c + bxx * x**2 + byy * y**2 + bzz * z**2 = 0
c C = cc + cx * x + cy * y + cz * z
c + cxy * x * y + cyz * y * z + czx * z * x
c + cxx * x**2 + cyy * y**2 + czz * z**2 = 0
c
c The three quadric surfaces have the normal vectors:
c
c AN = (ANx, ANy, ANz) = grad (A)
c ANx = ax + 2 * axx * x + axy * y + azx * z
c ANy = ay + axy * x + 2 * ayy * y + ayz * z
c ANz = az + azx * x + ayz * y + 2 * azz * z
c
c BN = (BNx, BNy, BNz) = grad (B)
c BNx = bx + 2 * bxx * x + bxy * y + bzx * z
c BNy = by + bxy * x + 2 * byy * y + byz * z
c BNz = bz + bzx * x + byz * y + 2 * bzz * z
c
c CN = (CNx, CNy, CNz) = grad (C)
c CNx = cx + 2 * cxx * x + cxy * y + czx * z
c CNy = cy + cxy * x + 2 * cyy * y + cyz * z
c CNz = cz + czx * x + cyz * y + 2 * czz * z
c
c The method is iterative, starting with the initial guess
c (x, y, z) = (px, py, pz). Each step of the iteration, the
c estimated triple point is successively replaced by its proximal
c point on each of the three surfaces (which works well for
c surfaces that are highly curved or meet nearly orthogonally),
c then improved by solving the three simultaneous equations:
c
c A + ANx * (x' - x) + ANy * (y' - y) + ANz * (z' - z) = 0
c B + BNx * (x' - x) + BNy * (y' - y) + BNz * (z' - z) = 0
c C + CNx * (x' - x) + CNy * (y' - y) + CNz * (z' - z) = 0
c
c (which works well for surfaces that are flat or meet at a small
c angle). Iterating is continued to convergence (nerr = 0) or
c failure (nerr > 0).
c
c Convergence is determined by the precision limit tol, which
c should be in the range from 1.e-5 to 1.e-15 for 64-bit
c floating point arithmetic.
c
c Input: px, py, pz,
c ac, ax, ay, az, axy, ayz, azx, axx, ayy, azz,
c bc, bx, by, bz, bxy, byz, bzx, bxx, byy, bzz,
c cc, cx, cy, cz, cxy, cyz, czx, cxx, cyy, czz, tol.
c
c Output: tx, ty, tz, nerr.
c
c Calls: aptsolv, aptwhis
c
c Glossary:
c
c ac, ... Input Coefficients of the implicit equation for the first
c quadric surface:
c ac + ax * x + ay * y + az * z +
c axy * x * y + ayz * y * z + azx * z * x +
c axx * x**2 + ayy * y**2 + azz * z**2 = 0
c
c bc, ... Input Coefficients of the implicit equation for the second
c quadric surface:
c bc + bx * x + by * y + bz * z +
c bxy * x * y + byz * y * z + bzx * z * x +
c bxx * x**2 + byy * y**2 + bzz * z**2 = 0
c
c cc, ... Input Coefficients of the implicit equation for the third
c quadric surface:
c cc + cx * x + cy * y + cz * z +
c cxy * x * y + cyz * y * z + czx * z * x +
c cxx * x**2 + cyy * y**2 + czz * z**2 = 0
c
c nerr Output Indicates an input error, if not zero:
c 1 if there is no solution, such as when there is no
c triple point.
c 2 if the method fails to converge in 1000 iterations,
c which may indicate too small a value of tol.
c 3 if the the method diverges.
c
c px,py,pz Input The x, y, z coordinates of point "p", the initial guess
c for the triple point.
c
c tol Input Numerical tolerance limit. With 64-bit floating point
c numbers, 1.e-5 to 1.e-11 is recommended.
c Used both as an absolute and relative convergence
c criterion. Changes in triple point coordinates less
c than 10 * tol are treated as insignificant.
c
c tx,ty,tz Output Coordinates x, y and z at one of the triple point
c intersections (usually the one nearest point "p")
c of the three quadric surfaces, if nerr = 0.
c If no solution is found, point "t" is returned as
c (-1.e99, -1.e99, -1.e99).
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Local variables.
common /ltripq/ abc(3,3) ! Coefficients of equations (row, column).
common /ltripq/ cba(3,3) ! Inverse of matrix abc (row, column).
common /ltripq/ d(3) ! Right-hand side constants.
common /ltripq/ dm(3) ! Right-hand side constants check.
common /ltripq/ dr(3) ! Right-hand side constants residuals.
common /ltripq/ work(3,4) ! Working matrix.
common /ltripq/ xyz(3) ! Solution.
common /ltripq/ xyzm(3) ! Solution check.
common /ltripq/ xyzr(3) ! Solution residuals.
common /ltripq/ dist(4) ! Distance to a proximal point.
common /ltripq/ xx(4) ! Coordinate x of a proximal point.
common /ltripq/ yy(4) ! Coordinate x of a proximal point.
common /ltripq/ zz(4) ! Coordinate x of a proximal point.
c.... Initialize.
nerr = 0
tx = -1.e99
ty = -1.e99
tz = -1.e99
c=======================================================================********
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptripq finding a triple point of three',
cbug & ' quadric surfaces.',
cbug & ' tol=',1pe13.5)
cbug 9902 format ('Initial guess:' /
cbug & ' px, py, pz =',1p3e22.14 )
cbug 9903 format ('Quadric ac =',1pe22.14 /
cbug & ' ax,ay,az =',1p3e22.14 /
cbug & ' axy,ayz,azx=',1p3e22.14 /
cbug & ' axx,ayy,azz=',1p3e22.14 )
cbug 9904 format ('Quadric bc =',1pe22.14 /
cbug & ' bx,by,bz =',1p3e22.14 /
cbug & ' bxy,byz,bzx=',1p3e22.14 /
cbug & ' bxx,byy,bzz=',1p3e22.14 )
cbug 9905 format ('Quadric cc =',1pe22.14 /
cbug & ' cx,cy,cz =',1p3e22.14 /
cbug & ' cxy,cyz,czx=',1p3e22.14 /
cbug & ' cxx,cyy,czz=',1p3e22.14 )
cbug write ( 3, 9901) tol
cbug write ( 3, 9902) px, py, pz
cbug write ( 3, 9903) ac, ax, ay, az, axy, ayz, azx, axx, ayy, azz
cbug write ( 3, 9904) bc, bx, by, bz, bxy, byz, bzx, bxx, byy, bzz
cbug write ( 3, 9905) cc, cx, cy, cz, cxy, cyz, czx, cxx, cyy, czz
cbugc***DEBUG ends.
c.... Initialize for the iteration.
x = px
y = py
z = pz
xold3 = x
yold3 = y
zold3 = z
aside = -1.e99
bside = -1.e99
cside = -1.e99
xcyc1 = -1.e99
ycyc1 = -1.e99
zcyc1 = -1.e99
xcyc2 = -1.e99
ycyc2 = -1.e99
zcyc2 = -1.e99
xcyc3 = -1.e99
ycyc3 = -1.e99
zcyc3 = -1.e99
c=======================================================================********
c.... Iterate to find the solution.
do 380 nit = 1, 1000 ! Iterate to find solution.
cbugcbugc***DEBUG begins.
cbugcbug 9910 format (/ i3,1x,'x,y,z =',1p3e22.14 )
cbugcbug write ( 3, 9910) nit, x, y, z
cbugcbugc***DEBUG ends.
c.... Try three successive proximal points.
ncon = 0 ! Number of sequential unchanged points.
ncyc = 0 ! Number of cycling points.
c.... Move to proximal point on first surface.
xold1 = x
yold1 = y
zold1 = z
tolx = 1.e-15
call aptwhis (x, y, z,
& ac, ax, ay, az, axy, ayz, azx, axx, ayy, azz,
& tolx,
& dside, anx, any, anz,
& nprox, xx, yy, zz, dist, nerra)
if ((nerra .ne. 0) .or.
& (nprox .le. 0)) then
nerr = 1
go to 999
endif
nr = 1.0 + ranf () * nprox
x = xx(nr)
y = yy(nr)
z = zz(nr)
cbugcbugc***DEBUG begins.
cbugcbug 9911 format (i3,1x,'x,y,z 1 =',1p3e22.14 )
cbugcbug write ( 3, 9911) nit, x, y, z
cbugcbugc***DEBUG ends.
if ((abs (x - xold3) .le. tol * (10.0 + abs (xold3))) .and.
& (abs (y - yold3) .le. tol * (10.0 + abs (yold3))) .and.
& (abs (z - zold3) .le. tol * (10.0 + abs (zold3)))) then
ncon = ncon + 1
endif
if ((abs (x - xcyc1) .le. tol * abs (xcyc1)) .and.
& (abs (y - ycyc1) .le. tol * abs (ycyc1)) .and.
& (abs (z - zcyc1) .le. tol * abs (zcyc1))) then
ncyc = ncyc + 1
endif
xcyc1 = x
ycyc1 = y
zcyc1 = z
c.... Move to proximal point on second surface.
xold2 = x
yold2 = y
zold2 = z
call aptwhis (x, y, z,
& bc, bx, by, bz, bxy, byz, bzx, bxx, byy, bzz,
& tolx,
& dside, bnx, bny, bnz,
& nprox, xx, yy, zz, dist, nerra)
if ((nerra .ne. 0) .or.
& (nprox .le. 0)) then
nerr = 1
go to 999
endif
nr = 1.0 + ranf () * nprox
x = xx(nr)
y = yy(nr)
z = zz(nr)
cbugcbugc***DEBUG begins.
cbugcbug 9912 format (i3,1x,'x,y,z 2 =',1p3e22.14 )
cbugcbug write ( 3, 9912) nit, x, y, z
cbugcbugc***DEBUG ends.
if ((abs (x - xold1) .le. tol * (10.0 + abs (xold1))) .and.
& (abs (y - yold1) .le. tol * (10.0 + abs (yold1))) .and.
& (abs (z - zold1) .le. tol * (10.0 + abs (zold1)))) then
ncon = ncon + 1
endif
if ((abs (x - xcyc2) .le. tol * abs (xcyc2)) .and.
& (abs (y - ycyc2) .le. tol * abs (ycyc2)) .and.
& (abs (z - zcyc2) .le. tol * abs (zcyc2))) then
ncyc = ncyc + 1
endif
xcyc2 = x
ycyc2 = y
zcyc2 = z
c.... Move to proximal point on third surface.
xold3 = x
yold3 = y
zold3 = z
call aptwhis (x, y, z,
& cc, cx, cy, cz, cxy, cyz, czx, cxx, cyy, czz,
& tolx,
& dside, cnx, cny, cnz,
& nprox, xx, yy, zz, dist, nerra)
if ((nerra .ne. 0) .or.
& (nprox .le. 0)) then
nerr = 1
go to 999
endif
nr = 1.0 + ranf () * nprox
x = xx(nr)
y = yy(nr)
z = zz(nr)
cbugcbugc***DEBUG begins.
cbugcbug 9913 format (i3,1x,'x,y,z 3 =',1p3e22.14 )
cbugcbug write ( 3, 9913) nit, x, y, z
cbugcbugc***DEBUG ends.
if ((abs (x - xold2) .le. tol * (10.0 + abs (xold2))) .and.
& (abs (y - yold2) .le. tol * (10.0 + abs (yold2))) .and.
& (abs (z - zold2) .le. tol * (10.0 + abs (zold2)))) then
ncon = ncon + 1
endif
if ((abs (x - xcyc3) .le. tol * abs (xcyc3)) .and.
& (abs (y - ycyc3) .le. tol * abs (ycyc3)) .and.
& (abs (z - zcyc3) .le. tol * abs (zcyc3))) then
ncyc = ncyc + 1
endif
xcyc3 = x
ycyc3 = y
zcyc3 = z
if (ncon .eq. 3) then
cbugcbugc***DEBUG begins.
cbugcbug 9914 format ('Converged in',i3,' iterations.')
cbugcbug write ( 3, 9914) nit
cbugcbugc***DEBUG ends.
go to 410
endif
if (ncyc .eq. 3) then
cbugcbugc***DEBUG begins.
cbugcbug 9915 format ('Detected cycling in',i3,' iterations.')
cbugcbug write ( 3, 9915) nit
cbugcbugc***DEBUG ends.
nerr = 3
go to 999
endif
c.... Find the new values of the implicit equations.
asideold = aside
bsideold = bside
csideold = cside
aside = ac + ax * x + ay * y + az * z
& + axy * x * y + ayz * y * z + azx * z * x
& + axx * x**2 + ayy * y**2 + azz * z**2
bside = bc + bx * x + by * y + bz * z
& + bxy * x * y + byz * y * z + bzx * z * x
& + bxx * x**2 + byy * y**2 + bzz * z**2
cside = cc + cx * x + cy * y + cz * z
& + cxy * x * y + cyz * y * z + czx * z * x
& + cxx * x**2 + cyy * y**2 + czz * z**2
cbugcbugc***DEBUG begins.
cbugcbug 9921 format ('side(a,b,c) =',1p3e22.14)
cbugcbug write ( 3, 9921) aside, bside, cside
cbugcbugc***DEBUG ends.
c.... Truncate meaningless residuals to zero.
aserr = tol * (abs(ac) +
& 2*abs(ax*x) + 2*abs(ay*y) + 2*abs(az*z) +
& 3*abs(axy*x*y) + 3*abs(ayz*y*z) + 3*abs(azx*z*x) +
& 3*abs(axx*x**2) + 3*abs(ayy*y**2) + 3*abs(azz*z**2))
if (abs (aside) .le. aserr) aside = 0.0
bserr = tol * (abs(bc) +
& 2*abs(bx*x) + 2*abs(by*y) + 2*abs(bz*z) +
& 3*abs(bxy*x*y) + 3*abs(byz*y*z) + 3*abs(bzx*z*x) +
& 3*abs(axx*x**2) + 3*abs(ayy*y**2) + 3*abs(azz*z**2))
if (abs (bside) .le. bserr) bside = 0.0
cserr = tol * (abs(cc) +
& 2*abs(cx*x) + 2*abs(cy*y) + 2*abs(cz*z) +
& 3*abs(cxy*x*y) + 3*abs(cyz*y*z) + 3*abs(czx*z*x) +
& 3*abs(cxx*x**2) + 3*abs(cyy*y**2) + 3*abs(czz*z**2))
if (abs (cside) .le. cserr) cside = 0.0
cbugcbugc***DEBUG begins.
cbugcbug 9922 format ('serr(a,b,c) =',1p3e22.14)
cbugcbug 9923 format ('side(a,b,c)tr=',1p3e22.14)
cbugcbug write ( 3, 9922) aserr, bserr, cserr
cbugcbug write ( 3, 9923) aside, bside, cside
cbugcbugc***DEBUG ends.
c.... See if the test point is at a triple intersection point.
if ((aside .eq. 0.0) .and.
& (bside .eq. 0.0) .and.
& (cside .eq. 0.0)) then
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9914) nit
cbugcbugc***DEBUG ends.
go to 410
endif
c.... See if the method is not converging.
if ((abs (aside) .ge. abs (asideold)) .and.
& (abs (bside) .ge. abs (bsideold)) .and.
& (abs (cside) .ge. abs (csideold))) then
cbugcbugc***DEBUG begins.
cbugcbug 9931 format ('aptripq not converging.')
cbugcbug write ( 3, 9931)
cbugcbugc***DEBUG ends.
nerr = 3
go to 999
endif
c.... Find the components of the normal vectors.
anx = ax + 2 * axx * x + axy * y + azx * z
any = ay + axy * x + 2 * ayy * y + ayz * z
anz = az + azx * x + ayz * y + 2 * azz * z
bnx = bx + 2 * bxx * x + bxy * y + bzx * z
bny = by + bxy * x + 2 * byy * y + byz * z
bnz = bz + bzx * x + byz * y + 2 * bzz * z
cnx = cx + 2 * cxx * x + cxy * y + czx * z
cny = cy + cxy * x + 2 * cyy * y + cyz * z
cnz = cz + czx * x + cyz * y + 2 * czz * z
cbugcbugc***DEBUG begins.
cbugcbug 9941 format ('anx,y,z =',1p3e22.14)
cbugcbug 9942 format ('bnx,y,z =',1p3e22.14)
cbugcbug 9943 format ('cnx,y,z =',1p3e22.14)
cbugcbug write ( 3, 9941) anx, any, anz
cbugcbug write ( 3, 9942) bnx, bny, bnz
cbugcbug write ( 3, 9943) cnx, cny, cnz
cbugcbugc***DEBUG ends.
c.... Find the coefficients of the three equations to be solved.
abc(1,1) = anx
abc(2,1) = any
abc(3,1) = anz
d(1) = -aside
abc(1,2) = bnx
abc(2,2) = bny
abc(3,2) = bnz
d(2) = -bside
abc(1,3) = cnx
abc(2,3) = cny
abc(3,3) = cnz
d(3) = -cside
c.... Solve the three simultaneous equations.
call aptsolv (abc, d, 3, work, 3, 0., xyz, cba, xyzm, xyzr,
& dm, dr, nerra)
if (nerra .eq. 0) then
dx = xyz(1)
dy = xyz(2)
dz = xyz(3)
cbugcbugc***DEBUG begins.
cbugcbug 9951 format ('dx,dy,dz =',1p3e22.14)
cbugcbug write ( 3, 9951) dx, dy, dz
cbugcbugc***DEBUG ends.
x = x + dx
y = y + dy
z = z + dz
cbugcbugc***DEBUG begins.
cbugcbug 9955 format (i3,1x,'x,y,z new=',1p3e22.14 )
cbugcbug write ( 3, 9955) nit, x, y, z
cbugcbugc***DEBUG ends.
if ((abs (dx) .lt. tol * (10.0 + abs (x))) .and.
& (abs (dy) .lt. tol * (10.0 + abs (y))) .and.
& (abs (dz) .lt. tol * (10.0 + abs (z)))) then
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9914) nit
cbugcbugc***DEBUG ends.
go to 410
endif
else
cbugcbugc***DEBUG begins.
cbugcbug 9958 format ('aptsolv says singular equations.')
cbugcbug write ( 3, 9958)
cbugcbugc***DEBUG ends.
nerr = 1
go to 999
endif
380 continue ! End of each iteration.
cbugcbugc***DEBUG begins.
cbugcbug 9961 format ('aptripq failed to converge in',i4,' iterations.')
cbugcbug write ( 3, 9961) nit
cbugcbugc***DEBUG ends.
nerr = 2 ! Failed to converge.
go to 999
c.... Converged to triple point.
410 if (nerr .eq. 0) then
tx = x
ty = y
tz = z
endif
c=======================================================================********
999 continue
cbugc***DEBUG begins.
cbug 9991 format (/ 'aptripq results. nerr=',i2,', nit=',i3.'.')
cbug 9992 format ('Triple (xyz)=',1p3e22.14 )
cbug write ( 3, 9991) nerr, nit
cbug write ( 3, 9992) tx, ty, tz
cbugc***DEBUG ends.
return
c.... End of subroutine aptripq. (+1 line)
end
UCRL-WEB-209832