C Q(1) = S*H(1) Q(2) = S*H(2) Q(3) = S*H(3) P(1) = Q(2)*A(3) - Q(3)*A(2) P(2) = Q(3)*A(1) - Q(1)*A(3) P(3) = Q(1)*A(2) - Q(2)*A(1) QTG = Q(1)*G(1) + Q(2)*G(2) + Q(3)*G(3) C C R = R + S*[(A X G) + (H X F) + QTG*P] C R(1) = R(1) + S*(A(2)*G(3)-A(3)*G(2) + . H(2)*F(3)-H(3)*F(2) + QTG*P(1)) R(2) = R(2) + S*(A(3)*G(1)-A(1)*G(3) + . H(3)*F(1)-H(1)*F(3) + QTG*P(2)) R(3) = R(3) + S*(A(1)*G(2)-A(2)*G(1) + . H(1)*F(2)-H(2)*F(1) + QTG*P(3)) C C T = T + S*(*I - A*A**T - P*P**T) C A1S = A(1)*A(1) A2S = A(2)*A(2) A3S = A(3)*A(3) T(1) = T(1) + S*(A2S + A3S - P(1)*P(1)) T(2) = T(2) - S*(A(2)*A(1) + P(2)*P(1)) T(3) = T(3) - S*(A(3)*A(1) + P(3)*P(1)) T(4) = T(4) + S*(A1S + A3S - P(2)*P(2)) T(5) = T(5) - S*(A(3)*A(2) + P(3)*P(2)) T(6) = T(6) + S*(A1S + A2S - P(3)*P(3)) RETURN END SUBROUTINE SOR1 (KM,K,F,D1F,D2F,OMEGA, NIT,DMAX,SA1, . SA2, C,CTC,G,GTG,IER) INTEGER KM, K, NIT, IER DOUBLE PRECISION F(3,0:KM,0:KM), D1F(3,0:KM,0:KM), . D2F(3,0:KM,0:KM), OMEGA, DMAX, . SA1(KM,KM), SA2(KM,KM), . C(3,0:KM,0:KM), CTC, . G(3,0:KM,0:KM), GTG C C*********************************************************** C C From MINSURF C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C (817) 565-2816 C 04/16/94 C C This subroutine employs a block SOR (Successive Over- C relaxation) method (with 3 by 3 blocks) to solve the C symmetric positive definite linear system associated with C computing the F-gradient G of PHI (the Sobolev gradient): C C L(G) = .5*D1[( D2(F) X D1(G) X D2(F) + C (D2(G) X D1(F)) X D2(F) )/FN] + C .5*D2[( D1(F) X D2(G) X D1(F) + C (D1(G) X D2(F)) X D1(F) )/FN] = L(F) and C C G = 0 on the boundary of the unit square, C C where F and G are vectors of length 3, D1 and D2 are first C partial derivative operators, and FN is the Euclidean norm C of the surface normal D1(F) X D2(F). Note that C C L(F) = D1[( D2(F) X D1(F) X D2(F) )/FN] + C D2[( D1(F) X D2(F) X D1(F) )/FN] = C(F) C C is the negative L2-gradient of the surface area PHI(F), C and is proportional to the mean curvature vector of the C surface F. C C The initial solution estimate is taken to be G = 0. C C On input: C C KM = Value of K used in the dimension statements C in the calling program. C C K = Number of cells in each direction. The dis- C cretization is defined by partitioning the unit C square into K**2 square cells. 2 .LE. K .LE. C KM. C C F = Array dimensioned 3 by 0:KM by 0:KM containing C grid point values of the parametric surface F: C F(1,i,j), F(2,i,j), and F(3,i,j) are the com- C ponents of the surface at (i/K,j/K) for i,j = C 0 to K. C C D1F,D2F = Arrays dimensioned 3 by 0:KM by 0:KM con- C taining central difference approximations C to D1(F) and D2(F) at the edge centers C scaled by 1/K. Thus, omitting the first C subscript, C C D1F(i,j) = F(i,j) - F(i-1,j) C D2F(i,j) = F(i,j) - F(i,j-1) C C for i,j = 0 to K, where D1F(0,j) and C D2F(i,0) are not defined (or used). C C OMEGA = Relaxation factor for the SOR method. C OMEGA = 1 corresponds to the Gauss-Seidel C method. 1 .LE. OMEGA .LE. 2. C C The above parameters are not altered by this routine. C C NIT = Maximum number of SOR iterations to be em- C ployed. This maximum will likely be achieved C if DMAX is smaller than the machine precision. C NIT .GE. 0. If NIT = 0, the L2-gradient C -L(F) is returned in G. C C DMAX = Nonnegative convergence criterion. The C method is terminated when the maximum (over C interior grid point values i,j = 1 to K-1) C change in U = F-G between iterations is at C most DMAX. The change in U is taken to be C the max-norm of the difference relative to 1 C plus the max-norm of the old value. C C SA1,SA2 = Arrays dimensioned KM by KM containing C areas of surface triangles scaled by 2: C C SA1(i,j) = Abs(D1F(i,j) X D2F(i,j)) C SA2(i,j) = Abs(D1F(i,j-1) X D2F(i-1,j)) C C for i,j = 1 to K, where Abs denotes C Euclidean norm. These quantities must C all be positive. C C C,G = Arrays dimensioned 3 by 0:KM by 0:KM. C C On output: C C NIT = Number of SOR iterations employed unless C IER = -1. C C DMAX = Maximum relative change in a value of U at C the last iteration unless IER = -1 or IER = C -4. C C SA1,SA2 = Arrays overwritten with the reciprocals of C the triangle areas unless IER = -1 or C IER = -4 (arrays partially overwritten). C C C = Array containing grid-point values of the mean C curvature vector (negative L2-gradient) with C zeros on the boundary. C is not altered if C IER = -1 or IER = -4. C C CTC = Sum (over the grid points) of squares of the C Euclidean norms of the mean curvature vectors C unless IER = -1 or IER = -4. C C G = Array containing grid-point values of the F- C gradient of PHI(F) satisfying L(G) = L(F) or, C if NIT = 0 on input, the L2-gradient: G = C -L(F) = -C. In either case, G = 0 on the boun- C dary of the unit square -- G(1,i,j) = G(2,i,j) = C G(3,i,j) = 0 for i = 0, i = K, j = 0, or j = K. C G is not altered if IER = -1 or IER = -4. C C GTG = Sum (over the grid points) of squares of the C Euclidean norms of the grid-point F-gradient C (or L2-gradient) values unless IER < 0. C C IER = Error indicator: C IER = 0 if no errors were encountered and the C convergence criterion was achieved (or C NIT = 0 on input). C IER = 1 if no errors were encountered but con- C vergence was not achieved within NIT C iterations. C IER = -1 if KM, K, OMEGA, NIT, or DMAX is out- C side its valid range on input. C IER = -2 if the linear system is numerically C singular. C IER = -3 if GTG exceeds GTGMAX (an arbitrary C large number defined below) indica- C ting a nearly singular system. C IER = -4 if SA1(i,j) .LE. 0 or SA2(i,j) .LE. 0 C for some i,j in the range [1,K] X C [1,K]. C C MINSURF modules required by SOR1: MAT, RES1 C C Intrinsic functions called by SOR1: DABS, DMAX1 C C*********************************************************** C DOUBLE PRECISION A(6), A22, A32, A33, B(3), CSUM, . D(3), DET, DMX, DU(3), GSUM, GTGMAX, . OM, R(3), R2, R3, T(3), TOL, U(3) INTEGER I, IM1, IP1, ITER, ITMAX, J, JM1, JP1, KK, . KM1, L C DATA GTGMAX/1.D10/ C KK = K KM1 = KK - 1 OM = OMEGA ITMAX = NIT TOL = DMAX C C Test for error flag -1 and initialize the iteration count C ITER and sum of squares of mean curvatures CSUM. C IF (KK .LT. 2 .OR. KK .GT. KM .OR. OM .LT. 1.D0 . .OR. OM .GT. 2.D0 .OR. ITMAX .LT. 0 .OR. . TOL .LT. 0.D0) GO TO 31 ITER = 0 CSUM = 0.D0 C C Test for error flag -4 and overwrite SA1 and SA2 with C reciprocals. C DO 2 J = 1,KK DO 1 I = 1,KK IF (SA1(I,J) .LE. 0.D0 .OR. SA2(I,J) .LE. 0.D0) . GO TO 34 SA1(I,J) = 1.D0/SA1(I,J) SA2(I,J) = 1.D0/SA2(I,J) 1 CONTINUE 2 CONTINUE C C Copy F into G as the initial estimate of U = F-G, and C store zero boundary values in C. C DO 4 J = 0,KK DO 3 I = 0,KK G(1,I,J) = F(1,I,J) G(2,I,J) = F(2,I,J) G(3,I,J) = F(3,I,J) C(1,I,J) = 0.D0 C(2,I,J) = 0.D0 C(3,I,J) = 0.D0 3 CONTINUE 4 CONTINUE C C Top of iteration loop: initialize DMX, and loop on C interior grid points. C 5 DMX = 0.D0 DO 14 J = 1,KM1 JM1 = J-1 JP1 = J+1 DO 13 I = 1,KM1 IM1 = I-1 IP1 = I+1 C C Initialize components of the order-3 system for the C change DU in the ij-th solution components (symmetric C matrix in A and residual in R), and store U(i,j) in U. C DO 6 L = 1,3 A(L) = 0.D0 A(L+3) = 0.D0 R(L) = 0.D0 U(L) = G(L,I,J) 6 CONTINUE C C Residual component: R(i,j) = C C { D2F(i+1,j) X (U(i+1,j)-U(i,j)) X D2F(i+1,j) C -D2F(i+1,j) X [(U(i+1,j)-U(i+1,j-1)) X D1F(i+1,j)]}/ C SA1(i+1,j) C + { D2F(i,j+1) X (U(i+1,j)-U(i,j)) X D2F(i,j+1) C -D2F(i,j+1) X [(U(i,j+1)-U(i,j)) X D1F(i+1,j)] C +D1F(i+1,j) X (U(i,j+1)-U(i,j)) X D1F(i+1,j) C -D1F(i+1,j) X [(U(i+1,j)-U(i,j)) X D2F(i,j+1)]}/ C SA2(i+1,j+1) C - { D2F(i,j) X (U(i,j)-U(i-1,j)) X D2F(i,j) C -D2F(i,j) X [(U(i,j)-U(i,j-1)) X D1F(i,j)] C +D1F(i,j) X (U(i,j)-U(i,j-1)) X D1F(i,j) C -D1F(i,j) X [(U(i,j)-U(i-1,j)) X D2F(i,j)]}/ C SA1(i,j) C - { D2F(i-1,j+1) X (U(i,j)-U(i-1,j)) X D2F(i-1,j+1) C -D2F(i-1,j+1) X [(U(i-1,j+1)-U(i-1,j)) X D1F(i,j)]}/ C SA2(i,j+1) C + { D1F(i,j+1) X (U(i,j+1)-U(i,j)) X D1F(i,j+1) C -D1F(i,j+1) X [(U(i,j+1)-U(i-1,j+1)) X D2F(i,j+1)]}/ C SA1(i,j+1) C - { D1F(i+1,j-1) X (U(i,j)-U(i,j-1)) X D1F(i+1,j-1) C -D1F(i+1,j-1) X [(U(i+1,j-1)-U(i,j-1)) X D2F(i,j)]}/ C SA2(i+1,j) C DO 7 L = 1,3 B(L) = G(L,IP1,J) - U(L) D(L) = G(L,IP1,J) - G(L,IP1,JM1) 7 CONTINUE CALL MAT (D2F(1,IP1,J),SA1(IP1,J), A ) CALL RES1 (D2F(1,IP1,J),B,D2F(1,IP1,J),D, . D1F(1,IP1,J),SA1(IP1,J), R ) C DO 8 L = 1,3 T(L) = D2F(L,I,JP1) - D1F(L,IP1,J) D(L) = G(L,I,JP1) - U(L) 8 CONTINUE CALL MAT (T,SA2(IP1,JP1), A ) CALL RES1 (T,B,D2F(1,I,JP1),D, . D1F(1,IP1,J),SA2(IP1,JP1), R ) C DO 9 L = 1,3 T(L) = D2F(L,I,J) - D1F(L,I,J) B(L) = G(L,IM1,J) - U(L) D(L) = G(L,I,JM1) - U(L) 9 CONTINUE CALL MAT (T,SA1(I,J), A ) CALL RES1 (T,B,D2F(1,I,J),D, . D1F(1,I,J),SA1(I,J), R ) C DO 10 L = 1,3 D(L) = G(L,IM1,J) - G(L,IM1,JP1) 10 CONTINUE CALL MAT (D2F(1,IM1,JP1),SA2(I,JP1), A ) CALL RES1 (D2F(1,IM1,JP1),B,D2F(1,IM1,JP1),D, . D1F(1,I,J),SA2(I,JP1), R ) C DO 11 L = 1,3 B(L) = G(L,I,JP1) - U(L) D(L) = G(L,I,JP1) - G(L,IM1,JP1) 11 CONTINUE CALL MAT (D1F(1,I,JP1),SA1(I,JP1), A ) CALL RES1 (D1F(1,I,JP1),B,D1F(1,I,JP1),D, . D2F(1,I,JP1),SA1(I,JP1), R ) C DO 12 L = 1,3 B(L) = G(L,I,JM1) - U(L) D(L) = G(L,I,JM1) - G(L,IP1,JM1) 12 CONTINUE CALL MAT (D1F(1,IP1,JM1),SA2(IP1,J), A ) CALL RES1 (D1F(1,IP1,JM1),B,D1F(1,IP1,JM1),D, . D2F(1,I,J),SA2(IP1,J), R ) IF (ITER .EQ. 0) THEN C C Store C entries and update CSUM on the first iteration C (U = F). C C(1,I,J) = R(1) C(2,I,J) = R(2) C(3,I,J) = R(3) CSUM = CSUM + R(1)**2 + R(2)**2 + R(3)**2 ENDIF IF (ITMAX .EQ. 0) GO TO 13 C C Solve the order-3 system associated with block ij. The C lower triangle of the matrix is stored by columns in C A: (A11,A21,A31,A22,A32,A33). C A22 = A(1)*A(4) - A(2)*A(2) A32 = A(1)*A(5) - A(2)*A(3) A33 = A(1)*A(6) - A(3)*A(3) R2 = A(1)*R(2) - A(2)*R(1) R3 = A(1)*R(3) - A(3)*R(1) DET = A22*A33 - A32*A32 IF (DET .EQ. 0.D0 .OR. A22 .EQ. 0.D0 .OR. . A(1) .EQ. 0.D0) GO TO 32 DU(3) = (A22*R3 - A32*R2)/DET DU(2) = (R2 - A32*DU(3))/A22 DU(1) = (R(1) - A(2)*DU(2) - A(3)*DU(3))/A(1) C C Update the solution components U(i,j) and the maximum C relative change in U. C DET = DMAX1(DABS(U(1)),DABS(U(2)),DABS(U(3))) . + 1.D0 G(1,I,J) = U(1) + OM*DU(1) G(2,I,J) = U(2) + OM*DU(2) G(3,I,J) = U(3) + OM*DU(3) DMX = DMAX1(DMX,DMAX1(DABS(DU(1)),DABS(DU(2)), . DABS(DU(3)))/DET) 13 CONTINUE 14 CONTINUE IF (ITMAX .EQ. 0) THEN C C Return G = -C for ITMAX = 0. C DO 16 J = 0,KK DO 15 I = 0,KK G(1,I,J) = -C(1,I,J) G(2,I,J) = -C(2,I,J) G(3,I,J) = -C(3,I,J) 15 CONTINUE 16 CONTINUE NIT = 0 DMAX = 0.E0 CTC = CSUM GTG = CSUM IER = 0 ELSE C C Increment ITER and test for convergence. C ITER = ITER + 1 IF (DMX .GT. TOL .AND. ITER .LT. ITMAX) GO TO 5 C C Compute G = F-U and GSUM = GTG. C GSUM = 0.D0 DO 18 J = 0,KK DO 17 I = 0,KK G(1,I,J) = F(1,I,J) - G(1,I,J) G(2,I,J) = F(2,I,J) - G(2,I,J) G(3,I,J) = F(3,I,J) - G(3,I,J) GSUM = GSUM + G(1,I,J)**2 + G(2,I,J)**2 + . G(3,I,J)**2 17 CONTINUE IF (GSUM .GT. GTGMAX) GO TO 33 18 CONTINUE C NIT = ITER DMAX = DMX CTC = CSUM GTG = GSUM IF (DMX .LE. TOL) THEN IER = 0 ELSE IER = 1 ENDIF ENDIF RETURN C C Error -1: KM, K, OMEGA, NIT, or DMAX is outside its valid C range on input. C 31 IER = -1 RETURN C C Error -2: A singular block was encountered. C 32 NIT = ITER DMAX = DMX CTC = CSUM IER = -2 RETURN C C Error -3: GTG would have exceeded GTGMAX. C 33 NIT = ITER DMAX = DMX CTC = CSUM IER = -3 RETURN C C Error -4: SA1(i,j) .LE. 0 or SA2(i,j) .LE. 0. C 34 NIT = 0 IER = -4 RETURN END SUBROUTINE SOR2 (KM,K,F,D1F,D2F,OMEGA, NIT,DMAX,SA1, . SA2, C,CTC,G,GTG,IER) INTEGER KM, K, NIT, IER DOUBLE PRECISION F(3,0:KM,0:KM), D1F(3,0:KM,0:KM), . D2F(3,0:KM,0:KM), OMEGA, DMAX, . SA1(KM,KM), SA2(KM,KM), . C(3,0:KM,0:KM), CTC, . G(3,0:KM,0:KM), GTG C C*********************************************************** C C From MINSURF C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C (817) 565-2816 C 04/16/94 C C This subroutine employs a block SOR (Successive Over- C relaxation) method (with 3 by 3 blocks) to solve the C symmetric positive definite linear system associated with C computing the F-gradient G of PHI (the Hessian gradient): C C L(G) = D1[( D2(F) X D1(G) X D2(F) + C (D2(G) X D1(F)) X D2(F) + C (D2(F) X D1(F)) X D2(G) - C S*(D2(F) X D1(F) X D2(F)) )/FN] + C D2[( D1(F) X D2(G) X D1(F) + C (D1(G) X D2(F)) X D1(F) + C (D1(F) X D2(F)) X D1(G) - C S*(D1(F) X D2(F) X D1(F)) )/FN] = L(F) and C C G = 0 on the boundary of the unit square, C C where F and G are vectors of length 3, D1 and D2 are first C partial derivative operators, FN is the Euclidean norm of C the surface normal D1(F) X D2(F), and C C S =/FN**2. C C Note that C C L(F) = D1[( D2(F) X D1(F) X D2(F) )/FN] + C D2[( D1(F) X D2(F) X D1(F) )/FN] = C(F) C C is the negative L2-gradient of the surface area PHI(F), C and is proportional to the mean curvature vector of the C surface F. C C The initial solution estimate is taken to be G = 0. C C On input: C C KM = Value of K used in the dimension statements C in the calling program. C C K = Number of cells in each direction. The dis- C cretization is defined by partitioning the unit C square into K**2 square cells. 2 .LE. K .LE. C KM. C C F = Array dimensioned 3 by 0:KM by 0:KM containing C grid point values of the parametric surface F: C F(1,i,j), F(2,i,j), and F(3,i,j) are the com- C ponents of the surface at (i/K,j/K) for i,j = C 0 to K. C C D1F,D2F = Arrays dimensioned 3 by 0:KM by 0:KM con- C taining central difference approximations C to D1(F) and D2(F) at the edge centers C scaled by 1/K. Thus, omitting the first C subscript, C C D1F(i,j) = F(i,j) - F(i-1,j) C D2F(i,j) = F(i,j) - F(i,j-1) C C for i,j = 0 to K, where D1F(0,j) and C D2F(i,0) are not defined (or used). C C OMEGA = Relaxation factor for the SOR method. C OMEGA = 1 corresponds to the Gauss-Seidel C method. 1 .LE. OMEGA .LE. 2. C C The above parameters are not altered by this routine. C C NIT = Maximum number of SOR iterations to be em- C ployed. This maximum will likely be achieved C if DMAX is smaller than the machine precision. C NIT .GE. 0. If NIT = 0, the L2-gradient C -L(F) is returned in G. C C DMAX = Nonnegative convergence criterion. The C method is terminated when the maximum (over C interior grid point values i,j = 1 to K-1) C change in U = F-G between iterations is at C most DMAX. The change in U is taken to be C the max-norm of the difference relative to 1 C plus the max-norm of the old value. C C SA1,SA2 = Arrays dimensioned KM by KM containing C areas of surface triangles scaled by 2: C C SA1(i,j) = Abs(D1F(i,j) X D2F(i,j)) C SA2(i,j) = Abs(D1F(i,j-1) X D2F(i-1,j)) C C for i,j = 1 to K, where Abs denotes C Euclidean norm. These quantities must C all be positive. C C C,G = Arrays dimensioned 3 by 0:KM by 0:KM. C C On output: C C NIT = Number of SOR iterations employed unless C IER = -1. C C DMAX = Maximum relative change in a value of U at C the last iteration unless IER = -1 or IER = C -4. C C SA1,SA2 = Arrays overwritten with the reciprocals of C the triangle areas unless IER = -1 or C IER = -4 (arrays partially overwritten) or C IER = 2 (not altered). C C C = Array containing grid-point values of the mean C curvature vector (negative L2-gradient) with C zeros on the boundary. C is not altered if C IER = -1 or IER = -4. C C CTC = Sum (over the grid points) of squares of the C Euclidean norms of the mean curvature vectors C unless IER = -1 or IER = -4. C C G = Array containing grid-point values of the F- C gradient of PHI(F) satisfying L(G) = L(F) or, C if NIT = 0 on input, the L2-gradient: G = C -L(F) = -C. In either case, G = 0 on the boun- C dary of the unit square -- G(1,i,j) = G(2,i,j) = C G(3,i,j) = 0 for i = 0, i = K, j = 0, or j = K. C G is not altered if IER = -1 or IER = -4. C C GTG = Sum (over the grid points) of squares of the C Euclidean norms of the grid-point F-gradient C (or L2-gradient) values unless IER < 0 or C IER = 2. C C IER = Error indicator: C IER = 0 if no errors were encountered and the C convergence criterion was achieved (or C NIT = 0 on input). C IER = 1 if no errors were encountered but con- C vergence was not achieved within NIT C iterations. C IER = 2 if the residual norm increased between C a pair of iterations, possibly due to C a non-positive definite Hessian. C IER = -1 if KM, K, OMEGA, NIT, or DMAX is out- C side its valid range on input. C IER = -2 if the linear system is numerically C singular. C IER = -3 if GTG exceeds GTGMAX (an arbitrary C large number defined below) indica- C ting a nearly singular system. C IER = -4 if SA1(i,j) .LE. 0 or SA2(i,j) .LE. 0 C for some i,j in the range [1,K] X C [1,K]. C C MINSURF module required by SOR2: RES2 C C Intrinsic functions called by SOR2: DABS, DMAX1 C C*********************************************************** C DOUBLE PRECISION A(6), A22, A32, A33, B(3), CSUM, . D(3), DET, DMX, DU(3), E(3), GSUM, . GTGMAX, OM, R(3), R2, R3, RTR, RTR0, . T(3), TOL, U(3) INTEGER I, IM1, IP1, ITER, ITMAX, J, JM1, JP1, KK, . KM1, L C DATA GTGMAX/1.D10/ C KK = K KM1 = KK - 1 OM = OMEGA ITMAX = NIT TOL = DMAX C C Test for error flag -1 and initialize the iteration count C ITER and sum of squares of mean curvatures CSUM. C IF (KK .LT. 2 .OR. KK .GT. KM .OR. OM .LT. 1.D0 . .OR. OM .GT. 2.D0 .OR. ITMAX .LT. 0 .OR. . TOL .LT. 0.D0) GO TO 31 ITER = 0 CSUM = 0.D0 C C Test for error flag -4 and overwrite SA1 and SA2 with C reciprocals. C DO 2 J = 1,KK DO 1 I = 1,KK IF (SA1(I,J) .LE. 0.D0 .OR. SA2(I,J) .LE. 0.D0) . GO TO 34 SA1(I,J) = 1.D0/SA1(I,J) SA2(I,J) = 1.D0/SA2(I,J) 1 CONTINUE 2 CONTINUE C C Copy F into G as the initial estimate of U = F-G, C store zero boundary values in C, and initialize RTR. C DO 4 J = 0,KK DO 3 I = 0,KK G(1,I,J) = F(1,I,J) G(2,I,J) = F(2,I,J) G(3,I,J) = F(3,I,J) C(1,I,J) = 0.D0 C(2,I,J) = 0.D0 C(3,I,J) = 0.D0 3 CONTINUE 4 CONTINUE RTR = GTGMAX C C Top of iteration loop: update RTR0 to the previous value C of RTR, initialize RTR and DMX, C and loop on interior grid points. C 5 RTR0 = RTR RTR = 0.D0 DMX = 0.D0 DO 14 J = 1,KM1 JM1 = J-1 JP1 = J+1 DO 13 I = 1,KM1 IM1 = I-1 IP1 = I+1 C C Initialize components of the order-3 system for the C change DU in the ij-th solution components (symmetric C matrix in A and residual in R), and store U(i,j) in U. C DO 6 L = 1,3 A(L) = 0.D0 A(L+3) = 0.D0 R(L) = 0.D0 U(L) = G(L,I,J) 6 CONTINUE C C Residual component: C C R(i,j) = Sum {s*[a X (b X c - d X e) + (c X e) X f + C s*s* *(c X e) X a]} C C Coefficient of U(i,j) -- symmetric positive semi-definite C order-3 matrix: C C A = Sum {s*[I - a*a**T - p*p**T]}, p = s*(c X e) X a C C where the sums are over the six triangles containing grid C point (i,j), and the scalar s and vectors a, b, c, d, e, C and f are defined as follows: C C [s, a, b, c, d, e, f] = C C 1) [SA1(i+1,j), D2F(i+1,j), U(i+1,j)-U(i,j), C D2F(i+1,j), U(i+1,j)-U(i+1,j-1), D1F(i+1,j), C U(i+1,j)-U(i+1,j-1)] C C 2) [SA2(i+1,j+1), D2F(i,j+1)-D1F(i+1,j), U(i+1,j)-U(i,j), C D2F(i,j+1), U(i,j+1)-U(i,j), D1F(i+1,j), C U(i,j+1)-U(i+1,j)] C C 3) [SA1(i,j), D2F(i,j)-D1F(i,j), U(i-1,j)-U(i,j), C D2F(i,j), U(i,j-1)-U(i,j), D1F(i,j), C U(i,j-1)-U(i-1,j)] C C 4) [SA2(i,j+1), D2F(i-1,j+1), U(i-1,j)-U(i,j), C D2F(i-1,j+1), U(i-1,j)-U(i-1,j+1), D1F(i,j), C U(i-1,j)-U(i-1,j+1)] C C 5) [SA1(i,j+1), D1F(i,j+1), U(i,j+1)-U(i,j), C D1F(i,j+1), U(i,j+1)-U(i-1,j+1), D2F(i,j+1), C U(i,j+1)-U(i-1,j+1)] C C 6) [SA2(i+1,j), D1F(i+1,j-1), U(i,j-1)-U(i,j), C D1F(i+1,j-1), U(i,j-1)-U(i+1,j-1), D2F(i,j), C U(i,j-1)-U(i+1,j-1)] C DO 7 L = 1,3 B(L) = G(L,IP1,J) - U(L) D(L) = G(L,IP1,J) - G(L,IP1,JM1) 7 CONTINUE CALL RES2 (D2F(1,IP1,J),B,D2F(1,IP1,J),D, . D1F(1,IP1,J),D,SA1(IP1,J), R,A ) C DO 8 L = 1,3 T(L) = D2F(L,I,JP1) - D1F(L,IP1,J) D(L) = G(L,I,JP1) - U(L) E(L) = G(L,I,JP1) - G(L,IP1,J) 8 CONTINUE CALL RES2 (T,B,D2F(1,I,JP1),D, . D1F(1,IP1,J),E,SA2(IP1,JP1), R,A ) C DO 9 L = 1,3 T(L) = D2F(L,I,J) - D1F(L,I,J) B(L) = G(L,IM1,J) - U(L) D(L) = G(L,I,JM1) - U(L) E(L) = G(L,I,JM1) - G(L,IM1,J) 9 CONTINUE CALL RES2 (T,B,D2F(1,I,J),D, . D1F(1,I,J),E,SA1(I,J), R,A ) C DO 10 L = 1,3 D(L) = G(L,IM1,J) - G(L,IM1,JP1) 10 CONTINUE CALL RES2 (D2F(1,IM1,JP1),B,D2F(1,IM1,JP1),D, . D1F(1,I,J),D,SA2(I,JP1), R,A ) C DO 11 L = 1,3 B(L) = G(L,I,JP1) - U(L) D(L) = G(L,I,JP1) - G(L,IM1,JP1) 11 CONTINUE CALL RES2 (D1F(1,I,JP1),B,D1F(1,I,JP1),D, . D2F(1,I,JP1),D,SA1(I,JP1), R,A ) C DO 12 L = 1,3 B(L) = G(L,I,JM1) - U(L) D(L) = G(L,I,JM1) - G(L,IP1,JM1) 12 CONTINUE CALL RES2 (D1F(1,IP1,JM1),B,D1F(1,IP1,JM1),D, . D2F(1,I,J),D,SA2(IP1,J), R,A ) IF (ITER .EQ. 0) THEN C C Store C entries and update CSUM on the first iteration C (U = F). C C(1,I,J) = R(1) C(2,I,J) = R(2) C(3,I,J) = R(3) CSUM = CSUM + R(1)**2 + R(2)**2 + R(3)**2 ENDIF IF (ITMAX .EQ. 0) GO TO 13 C C Update RTR. C RTR = RTR + R(1)**2 + R(2)**2 + R(3)**2 C C Solve the order-3 system associated with block ij. The C lower triangle of the matrix is stored by columns in C A: (A11,A21,A31,A22,A32,A33). C A22 = A(1)*A(4) - A(2)*A(2) A32 = A(1)*A(5) - A(2)*A(3) A33 = A(1)*A(6) - A(3)*A(3) R2 = A(1)*R(2) - A(2)*R(1) R3 = A(1)*R(3) - A(3)*R(1) DET = A22*A33 - A32*A32 IF (DET .EQ. 0.D0 .OR. A22 .EQ. 0.D0 .OR. . A(1) .EQ. 0.D0) GO TO 32 DU(3) = (A22*R3 - A32*R2)/DET DU(2) = (R2 - A32*DU(3))/A22 DU(1) = (R(1) - A(2)*DU(2) - A(3)*DU(3))/A(1) C C Update the solution components U(i,j) and the maximum C relative change in U. C DET = DMAX1(DABS(U(1)),DABS(U(2)),DABS(U(3))) . + 1.D0 G(1,I,J) = U(1) + OM*DU(1) G(2,I,J) = U(2) + OM*DU(2) G(3,I,J) = U(3) + OM*DU(3) DMX = DMAX1(DMX,DMAX1(DABS(DU(1)),DABS(DU(2)), . DABS(DU(3)))/DET) 13 CONTINUE 14 CONTINUE IF (ITMAX .EQ. 0) THEN C C Return G = -C for ITMAX = 0. C DO 16 J = 0,KK DO 15 I = 0,KK G(1,I,J) = -C(1,I,J) G(2,I,J) = -C(2,I,J) G(3,I,J) = -C(3,I,J) 15 CONTINUE 16 CONTINUE NIT = 0 DMAX = 0.E0 CTC = CSUM GTG = CSUM IER = 0 ELSE C C Increment ITER, test for an increase in RTR, and test for C convergence. C ITER = ITER + 1 IF (RTR .GT. RTR0) GO TO 19 IF (DMX .GT. TOL .AND. ITER .LT. ITMAX) GO TO 5 C C Compute G = F-U and GSUM = GTG. C GSUM = 0.D0 DO 18 J = 0,KK DO 17 I = 0,KK G(1,I,J) = F(1,I,J) - G(1,I,J) G(2,I,J) = F(2,I,J) - G(2,I,J) G(3,I,J) = F(3,I,J) - G(3,I,J) GSUM = GSUM + G(1,I,J)**2 + G(2,I,J)**2 + . G(3,I,J)**2 17 CONTINUE IF (GSUM .GT. GTGMAX) GO TO 33 18 CONTINUE C C Store the remaining output parameters. C NIT = ITER DMAX = DMX CTC = CSUM GTG = GSUM IF (DMX .LE. TOL) THEN IER = 0 ELSE IER = 1 ENDIF ENDIF RETURN C C Error 2: increase in RTR at iteration ITER. C 19 NIT = ITER DMAX = DMX CTC = CSUM C C Restore SA1 and SA2. C DO 21 J = 1,KK DO 20 I = 1,KK SA1(I,J) = 1.D0/SA1(I,J) SA2(I,J) = 1.D0/SA2(I,J) 20 CONTINUE 21 CONTINUE IER = 2 RETURN C C Error -1: KM, K, OMEGA, NIT, or DMAX is outside its valid C range on input. C 31 IER = -1 RETURN C C Error -2: a singular block was encountered. C 32 NIT = ITER DMAX = DMX CTC = CSUM IER = -2 RETURN C C Error -3: GTG would have exceeded GTGMAX. C 33 NIT = ITER DMAX = DMX CTC = CSUM IER = -3 RETURN C C Error -4: SA1(i,j) .LE. 0 or SA2(i,j) .LE. 0. C 34 NIT = 0 IER = -4 RETURN END