C C MINSURF (Parametric Minimal Surface Package) C 04/16/94 C C C Given a parametric surface FB(X,Y) = (F1,F2,F3) defined C on [0,1] X [0,1], this package finds a surface F that C agrees with FB on the boundary of the unit square and has C minimum surface area. Subroutine FCN is called to evalu- C ate FB (which provides an initial solution estimate as C well as the boundary values). C C The Fletcher-Reeves conjugate gradient method is used to C minimize the surface area functional C C PHI(F) = I( Abs(D1(F) X D2(F)) ) , C C where the cross product of first partials is the surface C normal vector, Abs denotes Euclidean norm, and I denotes C the integral over the unit square. Refer to Function PHI. C In place of the standard L2-gradient of PHI (which is C proportional to the negative mean curvature vector -C(F)), C the method employs an F-gradient G(F) corresponding to an C inner product that depends on the evaluation point F. It C is thus a variable metric method. Refer to Subroutine C SOR1 for a description of the Sobolev gradient and SOR2 C for the Hessian gradient. C C The unit square is partitioned into K**2 square cells C and each cell is partitioned into a pair of triangles by C the diagonal with slope -1. A triangle-based piecewise C linear finite element method is then employed to compute C values of (the 3-vector) F at the (K+1)**2 grid points. C C The following parameters may be altered by changing the C PARAMETER and DATA statements below. (If INTER = TRUE, C some parameters may be altered at run time.) C C CTOL = Nonnegative tolerance which, along with PHTOL, C defines convergence of the descent method if INTER C = FALSE: bound on the root-mean-square L2-gradient C of PHI. If CTOL < EPS, the tolerance is taken to C be EPS, where EPS = 100*(Machine precision). C C H0 = Default constant step-size if SEARCH = FALSE, or C initial value for the line search otherwise. C C INTER = Flag with value TRUE if and only if the program is C to be run interactively, in which case the user is C prompted for values of K, OMEGA, MAXCG, SEARCH (if C MAXCG = 0), and H0 (if MAXCG = 0 and SEARCH = C FALSE). C C K = Number of cells in each direction defining the dis- C cretization. The discretization error is proportional C to 1/K**2. K is stored by the first executable state- C ment below. C C KM = Maximum allowable value of K. The array storage C requirement in bytes is 184*KM**2 + 336*KM + 168 = C 476,968 for KM = 50. Note that a change in KM must C also be made in Subprograms LNSRCH and PHVAL. C C MAXCG = Number of Fletcher-Reeves conjugate gradient steps C to be taken before restarting with a steepest C descent step. C C NFMAX = Number of F-gradient steps to be taken before C switching to the standard L2-gradient if INTER = C FALSE. 1 .LE. NFMAX .LE. NGMAX. C C NGMAX = Maximum number of descent steps (outer iterations) C before termination without convergence if INTER = C FALSE. NGMAX .GE. 1. C C NPRNT = Number of descent steps between print steps. For C each print step, several lines of statistics de- C scribing the step are written to unit LPRT. The C first and last descent steps are necessarily print C steps, and NPRNT is overridden if INTER = TRUE. C NPRNT .GE. 1 C C NITMIN = Minimum number of SOR iterations for which SORTOL C is not decreased. Refer to STOL0 below. C C NITSOR = Maximum number of SOR iterations for each evalua- C tion of the F-gradient G. C C OMEGA = Relaxation factor in the range [1,2] for the SOR C method. If OMEGA = 0, a reasonable value is com- C puted at run time. C C PHTOL = Tolerance which, along with CTOL, defines conver- C gence of the descent method if INTER = FALSE: C bound on the relative change in PHI between iter- C ations. If either criterion is satisfied, the C method is terminated. If PHTOL < 0, convergence C is defined by CTOL. C C SEARCH = Flag with value TRUE if and only if a line search C for the optimal step-size is to be used at each C step of the descent method with the Sobolev grad- C ient. If SEARCH = FALSE, the step-size is taken C to be H0. C C STOL0 = Initial value of the nonnegative tolerance SORTOL C defining convergence of the SOR method: bound on C the maximum relative change in a solution compo- C nent between iterations. SORTOL is decreased by C a factor of 10 (but bounded below by CTOL) after C each call to SOR1 or SOR2 in which the number of C iterations is less than NITMIN. C C WRITF0 = Flag with value TRUE if and only if the initial C solution estimate F0 is to be written to disk. C Refer to LDSK below. C INTEGER KM PARAMETER (KM = 50) CHARACTER*80 TEXT DOUBLE PRECISION DSTORE, PHI DOUBLE PRECISION C(3,0:KM,0:KM), CTC, CTOL, CTOL2, . D(3,0:KM,0:KM), D1F(3,0:KM,0:KM), . D2F(3,0:KM,0:KM), DMAX, DPHI, EPS, . EPSP1, F(3,0:KM,0:KM), . F0(3,0:KM,0:KM), FMTOL, . G(3,0:KM,0:KM), GTGNEW, GTGOLD, H, . H0, OMEGA, PHI0, PHIF, PHTOL, Q, . RMSC, RMSG, RN, SA1(KM,KM), . SA2(KM,KM), SAMIN, SORTOL, STOL0, . X, Y INTEGER I, IER, IOPT, J, K, L, LDSK, LPRT, MAXCG, NCG, . NFMAX, NG, NGMAX, NIT, NITMIN, NITSOR, NITSM1, . NITSM2, NITSM3, NNEWT, NPRNT, NS LOGICAL INTER, NEWTON, PRNT, SEARCH, WRITF0 C C Common blocks for use by LNSRCH and PHVAL. C COMMON/PHCOM1/ D, K COMMON/PHCOM2/ F0 COMMON/PHCOM3/ F COMMON/PHCOM4/ D1F COMMON/PHCOM5/ D2F COMMON/PHCOM6/ SA1, SA2, SAMIN C DATA CTOL/1.D-13/, H0/1.D0/, INTER/.TRUE./, . MAXCG/2/, NFMAX/1000/, NGMAX/1000/, . NITMIN/10/, NITSOR/500/, NPRNT/1/, . OMEGA/0.D0/, PHTOL/.5D-4/, SEARCH/.TRUE./, . STOL0/1.D-3/, WRITF0/.FALSE./ C C LDSK = Logical unit number for writing the initial esti- C mate F0 (optional) and the solution F (in storage C order, three entries per line, with format 3D22.15: C left-to-right within bottom-to-top relative to the C grid). F0 and F are preceded by the value of K+1 C (format I4) on the first line. C C LPRT = Logical unit number for output other than the C solution. If INTER = TRUE, output is also written C to the screen. C DATA LDSK/2/, LPRT/3/ C C File names: C OPEN (LDSK,FILE='MINSURF.OUT') OPEN (LPRT,FILE='MINSURF.PRT') C C Interactive input formats: C 500 FORMAT (I5) 510 FORMAT (D22.15) 520 FORMAT (A80) C C Solution output formats: C 600 FORMAT (I4) 610 FORMAT (3D22.15) C C Noninteractive mode value of K: C K = 10 C C Compute a tolerance CTOL2 (bound on the root-mean-square C L2-gradient RMSC) defining the criterion for use of the C Hessian gradient (Newton's method) if INTER = FALSE. If C F is not sufficiently close to the solution, the Hessian C will not be positive definite, and the SOR method will C fail to converge after some number of wasted iterations. C In that case, the Sobolev gradient will be used. If C CTOL2 < CTOL, then no Newton steps will be used. This C parameter is not used if INTER = TRUE. C C The following expression (chosen on the basis of test C data) evaluates to 1.D-2 for K=10, 1.D-5 for K=30, and C 1.D-8 for K=50. C CTOL2 = 10.D0**(-.15D0*DBLE(K)-.5D0) C C Store a tolerance FMTOL for the line search (Subroutine C FMIN). FMTOL**2 should not be smaller than the machine C precision since the machine precision computed by FMIN C is incorrect on some systems. C C FMTOL = Sqrt(Machine Precision) C EPS = 1.D0 1 EPS = EPS/2.D0 EPSP1 = DSTORE(EPS + 1.D0) IF (EPSP1 .GT. 1.D0) GO TO 1 EPS = EPS*2.D0 FMTOL = DSQRT(EPS) C C CTOL is be bounded below by EPS = 100*(Machine Precision). C EPS = EPS*100.D0 IF (CTOL .LT. EPS) CTOL = EPS C IF (INTER) THEN C C Interactive mode: initialize NGMAX to 1. C NGMAX = 1 C C Print a message and prompt for a line of text to be C written to the print file. C WRITE (*,100) 100 FORMAT (///23X,'MINSURF: Minimal Surface Package'// . 10X,'Respond to each of the following ', . 'prompts with one of:'// . 15X,'a) the requested value followed by ', . 'Carriage Return,'/ . 15X,'b) Carriage Return only to specify ', . 'a default value, or'/ . 15X,'c) an invalid character (such as a ', . 'letter when numeric'/ . 19X,'input is requested) to back up to the', . ' previous option.'/// . 10X,'Specify a line of text with at most ', . '80 characters (the'/ . 10X,'current date for example) to be ', . 'included in MINSURF.PRT:'/) READ (*,520) TEXT C C Prompt the user for values of K, OMEGA, MAXCG, and WRITF0. C 2 WRITE (*,120) KM 120 FORMAT (//10X,'Specify the mesh size K in the ', . 'range 2 to ',I4/) READ (*,500,ERR=2) K IF (K .LT. 2 .OR. K .GT. KM) GO TO 2 C 3 WRITE (*,130) 130 FORMAT (//10X,'Specify the relaxation factor OMEGA', . ' in'/10X,'[1.,2.], or 0 to select the ', . 'default'/) READ (*,510,ERR=2) OMEGA IF (OMEGA .NE. 0.D0 .AND. (OMEGA .LT. 1.D0 .OR. . OMEGA .GT. 2.D0)) GO TO 3 C 4 SEARCH = .TRUE. H0 = 1.D0 WRITE (*,140) 140 FORMAT (//10X,'Specify the number of CG descent ', . 'steps MAXCG between'/ . 10X,'restarts with a steepest descent ', . 'iteration (MAXCG = 2'/ . 10X,'is recommended)'/) READ (*,500,ERR=3) MAXCG IF (MAXCG .LT. 0) GO TO 4 IF (MAXCG .EQ. 0) THEN C C Steepest Descent: prompt the user for the line search C option SEARCH. C WRITE (*,145) 145 FORMAT (//10X,'Use constant step-size (no line ', . 'search)? (0 ==> No, 1 ==> Yes)'/) READ (*,500,ERR=4) IOPT IF (IOPT .LT. 0 .OR. IOPT .GT. 1) GO TO 4 SEARCH = IOPT .EQ. 0 IF (.NOT. SEARCH) THEN C C No line search: prompt the user for the (constant) C step-size H0. C WRITE (*,148) H0 148 FORMAT (//10X,'Specify the step-size H0 > 0,'/ . 10X,'or 0 to select the default: ', . D7.1/) READ (*,510,ERR=4) X IF (X .LT. 0.D0) GO TO 4 IF (X .GT. 0.D0) H0 = X ENDIF ENDIF C 5 WRITE (*,150) 150 FORMAT (//10X,'Write the initial solution estimate ', . 'to disk?'/10X,'(0 ==> No, 1 ==> Yes)'/) READ (*,500,ERR=4) IOPT IF (IOPT .LT. 0 .OR. IOPT .GT. 1) GO TO 5 WRITF0 = IOPT .EQ. 1 ENDIF IF (OMEGA .EQ. 0.D0) THEN C C Compute a default value of OMEGA (optimal for Laplace's C equation). C X = .2D0*DBLE(K)**2 OMEGA = 2.D0*(1.D0+X)/(1.D0+X+DSQRT(1.D0+2.D0*X)) ENDIF C C Initialize SORTOL, and print a heading with parameter C values. C SORTOL = STOL0 WRITE (LPRT,160) 160 FORMAT (///27X,'MINSURF Output'/) IF (INTER) WRITE (LPRT,162) TEXT 162 FORMAT (5X,A80) WRITE (LPRT,164) K, SORTOL, NITSOR, OMEGA, MAXCG 164 FORMAT (//5X,'Mesh size (K**2 square cells): K = ',I4 . //5X,'SOR parameters:',2X, . 'Initial tolerance STOL0 = ',D7.1/ . 22X,'Max # iterations NITSOR = ',I5/ . 22X,'Relaxation factor OMEGA = ',F4.2// . 5X,'Number of CG iterations between ', . 'restarts: MAXCG = ',I3) IF (INTER) THEN WRITE (*,160) WRITE (*,162) TEXT WRITE (*,164) K, SORTOL, NITSOR, OMEGA, MAXCG ELSE WRITE (LPRT,170) CTOL, PHTOL, CTOL2, NGMAX, NFMAX 170 FORMAT (/5X,'Descent method parameters: ', . 'Convergence tol CTOL = ',D8.2/ . 31X,'Convergence tol PHTOL = ',D9.2/ . 29X,'Newton tolerance CTOL2 = ',D10.4/ . 35X,'Max # iterations NGMAX = ',I4/ . 23X,'Max # F-gradient evaluations NFMAX = ', . I4) ENDIF C C Initialize logical variables and iteration counts: C C NEWTON = TRUE iff the Hessian gradient is to be used C (rather than the Sobolev gradient), in which C case steepest descent iterations with step-size C 1 are Newton iterations. C C PRNT = TRUE iff statistics describing the iteration are C to be printed. C C NG = Number of outer iterations (gradient evaluations). C C NNEWT = Number of Newton iterations (included in NG). C C NITSM1 = Total number of SOR1 iterations. C C NITSM2 = Total number of SOR2 iterations. C C NITSM3 = Total number of wasted SOR2 iterations (in C calls to SOR2 without convergence). C NEWTON = .FALSE. PRNT = .TRUE. NG = 0 NNEWT = 0 NITSM1 = 0 NITSM2 = 0 NITSM3 = 0 C C Initialize F to F0, write F0 to disk (optionally), and C compute PHIF = PHI(F). C H = 1.D0/DBLE(K) DO 7 J = 0,K Y = DBLE(J)*H DO 6 I = 0,K X = DBLE(I)*H CALL FCN (X,Y, F(1,I,J)) 6 CONTINUE 7 CONTINUE IF (WRITF0) THEN WRITE (LDSK,600) K+1 WRITE (LDSK,610) (((F(L,I,J), L = 1,3),I = 0,K), . J = 0,K) ENDIF PHIF = PHI (KM,K,F, D1F,D2F,SA1,SA2,SAMIN) IF (SAMIN .LT. EPS) GO TO 21 C C Compute mean curvature vector C and Sobolev F-gradient C G at F. C C GTGOLD = G**T*G. C NIT = NITSOR DMAX = SORTOL CALL SOR1 (KM,K,F,D1F,D2F,OMEGA, NIT,DMAX,SA1,SA2, C, . CTC,G,GTGOLD,IER) NG = 1 NITSM1 = NIT IF (IER .LT. 0) GO TO 22 C C Compute root-mean-squares of C and G. C RN = DBLE((K+1)**2) RMSC = DSQRT(CTC/RN) RMSG = DSQRT(GTGOLD/RN) C C Print parameter values. C WRITE (LPRT,180) NG WRITE (LPRT,184) PHIF, NIT, DMAX, RMSC, RMSG IF (INTER) THEN WRITE (*,180) NG WRITE (*,184) PHIF, NIT, DMAX, RMSC, RMSG ENDIF 180 FORMAT (//15X,'Iteration ',I4,5X,'(Steepest Descent)') 181 FORMAT (//15X,'Iteration ',I4,5X,'(Newton step)') 182 FORMAT (//15X,'Iteration ',I4) 184 FORMAT (5X,'Surface area of F:',10X,'PHI(F) = ', . D21.15/ . 5X,'Number of SOR iterations: NIT = ', . I5/ . 5X,'Maximum change in U=F-G: DUMAX = ', . D8.2/ . 5X,'Root-mean-square curvature: RMSC = ', . D8.2/ . 5X,'Root-mean-square F-gradient: RMSG = ', . D8.2) C C Top of loop: start with a steepest descent iteration. C Initialize the number of CG iterations NCG, C set D = -G, and store F in F0. C 8 NCG = 0 DO 10 J = 0,K DO 9 I = 0,K D(1,I,J) = -G(1,I,J) D(2,I,J) = -G(2,I,J) D(3,I,J) = -G(3,I,J) F0(1,I,J) = F(1,I,J) F0(2,I,J) = F(2,I,J) F0(3,I,J) = F(3,I,J) 9 CONTINUE 10 CONTINUE C C Update F to F+H*D for optimal step-size H (or H = H0 if C SEARCH = FALSE or H = 1 if NEWTON = TRUE), and compute C PHI(F) and the relative change DPHI in PHI. C 11 PHI0 = PHIF IF (.NOT. NEWTON) THEN H = H0 CALL LNSRCH (PHI0,FMTOL,SEARCH, H, PHIF,IER) ELSE H = 1.D0 CALL LNSRCH (PHI0,FMTOL,.FALSE., H, PHIF,IER) ENDIF DPHI = (PHI0-PHIF)/PHI0 IF (PRNT) THEN C C Print SAMIN, H, PHIF, and DPHI. C WRITE (LPRT,190) H, SAMIN, PHIF, DPHI IF (INTER) WRITE (*,190) H, SAMIN, PHIF, DPHI 190 FORMAT (/5X,'Step-size:',23X,'H = ',D9.3/ . 5X,'Minimum triangle area: SAMIN = ', . D9.3/ . 5X,'Surface area: PHI(F+H*D) = ', . D21.15/ . 5X,'(PHI(F)-PHI(F+H*D))/PHI(F): DPHI = ', . D9.2/ . 1X,'______________________________________', . '___________________________') ENDIF IF (SAMIN .LT. EPS) GO TO 21 C C Test for termination. C IF (.NOT. INTER) THEN C C Noninteractive mode. C IF (RMSC .LE. CTOL .OR. DABS(DPHI) .LE. PHTOL . .OR. NG .GE. NGMAX) GO TO 30 ELSE C C Interactive mode. C IF (NG .GE. NGMAX) THEN 12 WRITE (*,200) 200 FORMAT (/10X,'Take another NS descent steps'/ . 10X,'(terminate if NS < 0) for NS = '/) READ (*,500,ERR=12) NS IF (NS .LT. 0) GO TO 30 IF (NS .EQ. 0) NS = 1 NGMAX = NGMAX + NS 13 WRITE (*,210) 210 FORMAT (10X,'Use the L2-gradient (in place of ', . 'the F-gradient)?'/ . 10X,'(0 ==> No, 1 ==> Yes)'/) READ (*,500,ERR=13) IOPT IF (IOPT .LT. 0 .OR. IOPT .GT. 1) GO TO 13 IF (IOPT .EQ. 0) THEN NFMAX = NGMAX NPRNT = 1 WRITE (*,220) 220 FORMAT (10X,'Try a Newton step? (0 ==> No, ', . '1 ==> Yes)'/) READ (*,500,ERR=13) IOPT IF (IOPT .LT. 0 .OR. IOPT .GT. 1) GO TO 13 NEWTON = IOPT .EQ. 1 ELSE NFMAX = 0 NPRNT = NGMAX ENDIF ENDIF ENDIF C C Compute mean curvature vector C and F-gradient G at F. C IF (.NOT. INTER .AND. RMSC .LT. CTOL2) . NEWTON = .TRUE. IF (NG .LT. NFMAX) THEN NIT = NITSOR ELSE NIT = 0 ENDIF 14 DMAX = SORTOL IF (NEWTON) THEN CALL SOR2 (KM,K,F,D1F,D2F,OMEGA, NIT,DMAX,SA1, . SA2, C,CTC,G,GTGNEW,IER) IF (IER .EQ. 2) THEN NITSM3 = NITSM3 + NIT NEWTON = .FALSE. WRITE (LPRT,230) NG+1, NIT, DMAX IF (INTER) WRITE (*,230) NG+1, NIT, DMAX 230 FORMAT (//5X,'Newton iteration (',I4,') failed:', . ' NIT =',I4,', DUMAX = ',D8.2) NIT = NITSOR GO TO 14 ELSE NNEWT = NNEWT + 1 NITSM2 = NITSM2 + NIT ENDIF ELSE CALL SOR1 (KM,K,F,D1F,D2F,OMEGA, NIT,DMAX,SA1, . SA2, C,CTC,G,GTGNEW,IER) NITSM1 = NITSM1 + NIT ENDIF NG = NG + 1 IF (IER .LT. 0) GO TO 22 IF (NIT .LT. NITMIN) . SORTOL = DMAX1(SORTOL/10.D0,CTOL) C C Compute root-mean-squares of C and G. C RMSC = DSQRT(CTC/RN) RMSG = DSQRT(GTGNEW/RN) PRNT = NPRNT*(NG/NPRNT) .EQ. NG IF (PRNT) THEN C C Print parameter values. C IF (NEWTON) THEN WRITE (LPRT,181) NG IF (INTER) WRITE (*,181) NG ELSEIF (NCG .EQ. MAXCG) THEN WRITE (LPRT,180) NG IF (INTER) WRITE (*,180) NG ELSE WRITE (LPRT,182) NG IF (INTER) WRITE (*,182) NG ENDIF WRITE (LPRT,184) PHIF, NIT, DMAX, RMSC, RMSG IF (INTER) WRITE (*,184) PHIF, NIT, DMAX, RMSC, RMSG ENDIF C C Compute a weight Q = GTGNEW/GTGOLD and update GTGOLD. C Q = GTGNEW/GTGOLD GTGOLD = GTGNEW IF (NCG .LT. MAXCG .AND. .NOT. NEWTON) THEN C C Compute a new search direction D = -G + Q*D, and C store F in F0. C NCG = NCG + 1 DO 16 J = 0,K DO 15 I = 0,K D(1,I,J) = -G(1,I,J) + Q*D(1,I,J) D(2,I,J) = -G(2,I,J) + Q*D(2,I,J) D(3,I,J) = -G(3,I,J) + Q*D(3,I,J) F0(1,I,J) = F(1,I,J) F0(2,I,J) = F(2,I,J) F0(3,I,J) = F(3,I,J) 15 CONTINUE 16 CONTINUE GO TO 11 ENDIF C C Bottom of loop. C GO TO 8 C C Error encountered in Function PHI. C 21 WRITE (LPRT,410) EPS IF (INTER) WRITE (*,410) EPS 410 FORMAT (///10X,'*** Error: triangle area less than ', . D9.3,' ***') GO TO 30 C C Error encountered in Subroutine SOR1 or SOR2. C 22 WRITE (LPRT,420) IER IF (INTER) WRITE (*,420) IER 420 FORMAT (///10X,'*** Error in SOR: IER = ',I2,' ***') C C Write the solution and terminate execution. C 30 WRITE (LDSK,600) K+1 WRITE (LDSK,610) (((F(L,I,J), L = 1,3), I = 0,K), . J = 0,K) IF (.NOT. PRNT) THEN C C Print parameter values. C IF (NEWTON) THEN WRITE (LPRT,181) NG IF (INTER) WRITE (*,181) NG ELSEIF (NCG .EQ. MAXCG) THEN WRITE (LPRT,180) NG IF (INTER) WRITE (*,180) NG ELSE WRITE (LPRT,182) NG IF (INTER) WRITE (*,182) NG ENDIF WRITE (LPRT,184) PHI0, NIT, DMAX, RMSC, RMSG IF (INTER) WRITE (*,184) PHI0, NIT, DMAX, RMSC, RMSG WRITE (LPRT,190) H, SAMIN, PHIF, DPHI IF (INTER) WRITE (*,190) H, SAMIN, PHIF, DPHI ENDIF C C Print iteration counts. C WRITE (LPRT,240) NG-NNEWT IF (INTER) WRITE (*,240) NG-NNEWT 240 FORMAT (//5X,'Number of CG or steepest descent ', . 'iterations: ',I4) WRITE (LPRT,245) NITSM1 IF (INTER) WRITE (*,245) NITSM1 245 FORMAT (/10X,'Total number of SOR1 iterations: ',I6) WRITE (LPRT,250) NNEWT IF (INTER) WRITE (*,250) NNEWT 250 FORMAT (/5X,'Number of Newton iterations: ',I4) WRITE (LPRT,255) NITSM2 IF (INTER) WRITE (*,255) NITSM2 255 FORMAT (/10X,'Total number of SOR2 iterations: ',I6) WRITE (LPRT,260) NITSM3 IF (INTER) WRITE (*,260) NITSM3 260 FORMAT (/10X,'Total number of unused SOR2 ', . 'iterations: ',I6) STOP END DOUBLE PRECISION FUNCTION DSTORE (X) DOUBLE PRECISION X C C*********************************************************** C C From TSPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C (817) 565-2767 C 06/01/90 C C This function forces its argument X to be stored in a C memory location, thus providing a means of determining C floating point number characteristics (such as the machine C precision) when it is necessary to avoid computation in C high precision registers. C C On input: C C X = Double precision value to be stored. C C X is not altered by this function. C C On output: C C STORE = Value of X after it has been stored and C possibly truncated or rounded to the double C precision word length. C C Modules required by STORE: None C C*********************************************************** C DOUBLE PRECISION Y COMMON/STCOM/Y Y = X DSTORE = Y RETURN END DOUBLE PRECISION FUNCTION FMIN(AX,BX,F,TOL) DOUBLE PRECISION AX,BX,F,TOL C***BEGIN PROLOGUE FMIN C***CATEGORY NO. G1A2 C***KEYWORD(S) ONE-DIMENSIONAL MINIMIZATION, UNIMODAL FUNCTION C***AUTHOR R. BRENT C***DATE WRITTEN 730101 (YYMMDD) C***PURPOSE C An approximation to the point where F attains a minimum on C the interval (AX,BX) is determined as the value of the function C FMIN. C C PARAMETERS C C INPUT C C AX (double precision) left endpoint of initial interval C BX (double precision) right endpoint of initial interval C F function subprogram which evaluates F(X) for any X C in the interval (AX,BX) C TOL (double precision) desired length of the interval of uncertainty C of the final result ( .ge. 0.0) C C C OUTPUT C C FMIN abcissa approximating the minimizer of F C AX lower bound for minimizer C BX upper bound for minimizer C C***DESCRIPTION C C The method used is a combination of golden section search and C successive parabolic interpolation. Convergence is never much C slower than that for a Fibonacci search. If F has a continuous C second derivative which is positive at the minimum (which is not C at AX or BX), then convergence is superlinear, and usually of the C order of about 1.324.... C C The function F is never evaluated at two points closer together C than EPS*ABS(FMIN) + (TOL/3), where EPS is approximately the C square root of the relative machine precision. If F is a unimodal C function and the computed values of F are always unimodal when C separated by at least EPS*ABS(XSTAR) + (TOL/3), then FMIN C approximates the abcissa of the global minimum of F on the C interval AX,BX with an error less than 3*EPS*ABS(FMIN) + TOL. C If F is not unimodal, then FMIN may approximate a local, but C perhaps non-global, minimum to the same accuracy. C C This function subprogram is a slightly modified version of the C ALGOL 60 procedure LOCALMIN given in Richard Brent, Algorithms for C Minimization Without Derivatives, Prentice-Hall, Inc. (1973). C C***REFERENCE(S) C Richard Brent, Algorithms for Minimization Without Derivatives, C Prentice-Hall, Inc. (1973). C***ROUTINES CALLED NONE C***END PROLOGUE DOUBLE PRECISION A,B,C,D,E,EPS,XM,P,Q,R,TOL1,TOL2,U,V,W DOUBLE PRECISION FU,FV,FW,FX,X DOUBLE PRECISION DABS,DSQRT,DSIGN EXTERNAL F C C C is the squared inverse of the golden ratio C C = 0.5D0*(3. - DSQRT(5.0D0)) C C EPS is approximately the square root of the relative machine C precision. C EPS = 1.0D00 10 EPS = EPS/2.0D00 TOL1 = 1.0D0 + EPS IF (TOL1 .GT. 1.0D00) GO TO 10 EPS = DSQRT(EPS) C C initialization C A = AX B = BX V = A + C*(B - A) W = V X = V E = 0.0D0 FX = F(X) FV = FX FW = FX C C main loop starts here C 20 XM = 0.5D0*(A + B) TOL1 = EPS*DABS(X) + TOL/3.0D0 TOL2 = 2.0D0*TOL1 C C check stopping criterion C IF (DABS(X - XM) .LE. (TOL2 - 0.5D0*(B - A))) GO TO 90 C C is golden-section necessary C IF (DABS(E) .LE. TOL1) GO TO 40 C C fit parabola C R = (X - W)*(FX - FV) Q = (X - V)*(FX - FW) P = (X - V)*Q - (X - W)*R Q = 2.0D00*(Q - R) IF (Q .GT. 0.0D0) P = -P Q = DABS(Q) R = E E = D C C is parabola acceptable C 30 IF (DABS(P) .GE. DABS(0.5D0*Q*R)) GO TO 40 IF (P .LE. Q*(A - X)) GO TO 40 IF (P .GE. Q*(B - X)) GO TO 40 C C a parabolic interpolation step C D = P/Q U = X + D C C F must not be evaluated too close to AX or BX C IF ((U - A) .LT. TOL2) D = DSIGN(TOL1, XM - X) IF ((B - U) .LT. TOL2) D = DSIGN(TOL1, XM - X) GO TO 50 C C a golden-section step C 40 IF (X .GE. XM) E = A - X IF (X .LT. XM) E = B - X D = C*E C C F must not be evaluated too close to X C 50 IF (DABS(D) .GE. TOL1) U = X + D IF (DABS(D) .LT. TOL1) U = X + DSIGN(TOL1, D) FU = F(U) C C update A, B, V, W, and X C IF (FU .GT. FX) GO TO 60 IF (U .GE. X) A = X IF (U .LT. X) B = X V = W FV = FW W = X FW = FX X = U FX = FU GO TO 20 60 IF (U .LT. X) A = U IF (U .GE. X) B = U IF (FU .LE. FW) GO TO 70 IF (W .EQ. X) GO TO 70 IF (FU .LE. FV) GO TO 80 IF (V .EQ. X) GO TO 80 IF (V .EQ. W) GO TO 80 GO TO 20 70 V = W FV = FW W = U FW = FU GO TO 20 80 V = U FV = FU GO TO 20 C C end of main loop C 90 FMIN = X RETURN END SUBROUTINE LNSRCH (PHI0,TOLFM,SEARCH, H, PHIF,IER) INTEGER IER LOGICAL SEARCH DOUBLE PRECISION PHI0, TOLFM, H, PHIF 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 03/28/94 C C This subroutine minimizes PHI(F0+H*D) over positive C step-sizes H, where PHI is the functional defined by Func- C tion PHI. C C On input: C C K = (PHCOM parameter) Number of cells in each direc- C tion. (There are (K+1)**2 grid points.) 2 .LE. C K .LE. KM. C C D = (PHCOM parameter) Array defining the search C direction for the descent method. C C F0 = (PHCOM parameter) Array containing values of C the current estimate of the minimizer of PHI. C C PHI0 = PHI(F0). C C TOLFM = Nonnegative tolerance for Subroutine FMIN: C the desired length of the interval of uncer- C tainty for the optimal step-size. C C SEARCH = Flag with value TRUE if and only if a line C search for the optimal step-size H is to be C employed. C C The above parameters are not altered by this subroutine. C C H = Initial estimate of the optimal step-size if C SEARCH = TRUE, or step-size to be used other- C wise. The value H = 1 or the output from a C previous call is a reasonable choice. H > 0. C C On output: C C H = Estimate of the optimal step-size if SEARCH = C TRUE, unaltered otherwise. C C F = (PHCOM parameter) Updated solution estimate: C F0+H*D. C C D1F,D2F,SA1,SA2,SAMIN = (PHCOM parameters) Output C values from a call to PHI C with F as input. C C PHIF = PHI(F). C C IER = Error indicator: C IER = 0 if no errors were encountered. C IER = 1 if K < 2, K > KM, TOLFM < 0, or H .LE. C 0 on input. Output parameters are not C altered in this case. C C Modules required by LNSRCH: FMIN, PHI, PHVAL C C Common blocks required by LNSRCH: PHCOM1, PHCOM2, PHCOM3, C PHCOM4, PHCOM5, PHCOM6 C C*********************************************************** C INTEGER K, KM PARAMETER (KM=50) DOUBLE PRECISION FMIN, PHVAL DOUBLE PRECISION D(3,0:KM,0:KM), F0(3,0:KM,0:KM), . F(3,0:KM,0:KM), D1F(3,0:KM,0:KM), . D2F(3,0:KM,0:KM), SA1(KM,KM), . SA2(KM,KM), SAMIN DOUBLE PRECISION A, B, PHB C C Common blocks used to pass parameters to PHVAL: C COMMON/PHCOM1/ D, K COMMON/PHCOM2/ F0 COMMON/PHCOM3/ F COMMON/PHCOM4/ D1F COMMON/PHCOM5/ D2F COMMON/PHCOM6/ SA1, SA2, SAMIN EXTERNAL PHVAL C C Test for invalid input. C IF (K .LT. 2 .OR. K .GT. KM .OR. TOLFM .LT. 0.D0 . .OR. H .LE. 0.D0) THEN IER = 1 RETURN ENDIF IF (SEARCH) THEN C C Find a bracketing interval [A,B] = [0,B] such that C PHI(F0+B*D) > PHI0 = PHI(F0). B is initialized to C 2*H and doubled at each step as necessary. C A = 0.D0 B = H 1 B = 2.D0*B PHB = PHVAL(B) IF (PHB .LE. PHI0) GO TO 1 C C Compute the optimal step-size H. C H = FMIN (A,B,PHVAL,TOLFM) ENDIF C C Compute the final value PHI(F0+H*D). C PHIF = PHVAL(H) IER = 0 RETURN END SUBROUTINE MAT (A,S, C ) DOUBLE PRECISION A(3), S, C(6) 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 01/08/92 C C Given a vector A of length 3, a scalar S, and a symmet- C ric matrix C, this subroutine overwrites C with C C C + S*(I - A*A**T) , C C where = (A**T)*A and A**T denotes the transpose of C A. The operation count is 9 adds and 12 multiples. C C If S > 0, the matrix added to C is symmetric and posi- C tive semi-definite and proportional to the orthogonal C projection onto the complement of A. A compact storage C scheme is used with C stored as a vector of length 6 C containing the lower triangle stored by columns, or equiv- C alently, the upper triangle stored by rows: C (C11, C21, C31, C22, C32, C33). C C Modules required by MAT: None C C*********************************************************** C DOUBLE PRECISION A1S, A2S, A3S C A1S = A(1)*A(1) A2S = A(2)*A(2) A3S = A(3)*A(3) C C(1) = C(1) + S*(A2S + A3S) C(2) = C(2) - S*A(2)*A(1) C(3) = C(3) - S*A(3)*A(1) C(4) = C(4) + S*(A1S + A3S) C(5) = C(5) - S*A(3)*A(2) C(6) = C(6) + S*(A1S + A2S) RETURN END DOUBLE PRECISION FUNCTION PHI (KM,K,F, D1F,D2F,SA1, . SA2,SAMIN) INTEGER KM, K DOUBLE PRECISION F(3,0:KM,0:KM), D1F(3,0:KM,0:KM), . D2F(3,0:KM,0:KM), SA1(KM,KM), . SA2(KM,KM), SAMIN 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/01/92 C C Given a K by K grid of square cells which partition the C unit square, and the grid point values (vectors of length C 3) of a parametric surface F, this function returns scaled C central difference approximations to the first partial C derivatives D1(F) and D2(F) at the edge centers, areas of C the triangular facets of a piecewise linear interpolant of C the grid point values, and the area PHI(F) of the piece- C wise linear interpolant on the triangulation defined by C partitioning each cell by the diagonal with slope -1. C C On input: C C KM = Value of K used in the dimension statements in C the calling program. C C K = Number of cells in each direction. (There are C (K+1)**2 grid points.) 2 .LE. K .LE. 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 On output: 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 altered). 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. C C SAMIN = Minimum of the 2*K**2 triangle areas. If C SAMIN = 0, F has a tangent plane discontin- C uity. C C PHI = Surface area of the piecewise linear interpo- C lant, or -1 if K < 2 or K > KM, in which case C the other output parameters are not defined. C C Modules required by PHI: None C C Intrinsic function called by PHI: DMIN1, DSQRT C C*********************************************************** C DOUBLE PRECISION FN1, FN2, FN3, SMIN, SUM INTEGER I, IM1, J, JM1, KK C C Test for invalid input. C KK = K IF (KK .LT. 2 .OR. KK .GT. KM) THEN PHI = -1.D0 RETURN ENDIF C C Store D1F(i,0) and D2F(0,j). C DO 1 I = 1,KK IM1 = I-1 D1F(1,I,0) = F(1,I,0) - F(1,IM1,0) D1F(2,I,0) = F(2,I,0) - F(2,IM1,0) D1F(3,I,0) = F(3,I,0) - F(3,IM1,0) D2F(1,0,I) = F(1,0,I) - F(1,0,IM1) D2F(2,0,I) = F(2,0,I) - F(2,0,IM1) D2F(3,0,I) = F(3,0,I) - F(3,0,IM1) 1 CONTINUE C C Initialize SMIN (smallest SA value encountered) and SUM C (sum of SA values). C SMIN = 1.D20 SUM = 0.D0 C C Store D1F(i,j), D2F(i,j), SA1(i,j), and SA2(i,j) for C i,j = 1 to K, and accumulate SMIN and SUM. C DO 3 J = 1,KK JM1 = J-1 DO 2 I = 1,KK IM1 = I-1 D1F(1,I,J) = F(1,I,J) - F(1,IM1,J) D1F(2,I,J) = F(2,I,J) - F(2,IM1,J) D1F(3,I,J) = F(3,I,J) - F(3,IM1,J) D2F(1,I,J) = F(1,I,J) - F(1,I,JM1) D2F(2,I,J) = F(2,I,J) - F(2,I,JM1) D2F(3,I,J) = F(3,I,J) - F(3,I,JM1) FN1 = D1F(2,I,J)*D2F(3,I,J)-D1F(3,I,J)*D2F(2,I,J) FN2 = D1F(3,I,J)*D2F(1,I,J)-D1F(1,I,J)*D2F(3,I,J) FN3 = D1F(1,I,J)*D2F(2,I,J)-D1F(2,I,J)*D2F(1,I,J) SA1(I,J) = DSQRT(FN1*FN1 + FN2*FN2 + FN3*FN3) FN1 = D1F(2,I,JM1)*D2F(3,IM1,J) - . D1F(3,I,JM1)*D2F(2,IM1,J) FN2 = D1F(3,I,JM1)*D2F(1,IM1,J) - . D1F(1,I,JM1)*D2F(3,IM1,J) FN3 = D1F(1,I,JM1)*D2F(2,IM1,J) - . D1F(2,I,JM1)*D2F(1,IM1,J) SA2(I,J) = DSQRT(FN1*FN1 + FN2*FN2 + FN3*FN3) SMIN = DMIN1(SMIN,SA1(I,J),SA2(I,J)) SUM = SUM + SA1(I,J) + SA2(I,J) 2 CONTINUE 3 CONTINUE C C Store SAMIN and PHI. C SAMIN = SMIN PHI = SUM/2.D0 RETURN END DOUBLE PRECISION FUNCTION PHVAL (H) DOUBLE PRECISION H 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 03/14/94 C C This function, designed to be passed as a parameter to C Function FMIN, returns the objective function value C PHI(F0+H*D) and updates common block parameters F = F0 + C H*D, D1F, D2F, SA1, SA2, and SAMIN. C C On input: C C H = Step-size in the direction D defining F. C C H is not altered by this function. C C On output: C C PHVAL = PHI(F), where F = F0+H*D. C C PHCOM parameters F, D1F, D2F, SA1, SA2, and SAMIN C are updated. C C Module required by PHVAL: PHI C C Common blocks required by PHVAL: PHCOM1, PHCOM2, PHCOM3, C PHCOM4, PHCOM5, PHCOM6 C C*********************************************************** C INTEGER I, J, K, KM PARAMETER (KM=50) DOUBLE PRECISION PHI DOUBLE PRECISION D(3,0:KM,0:KM), F0(3,0:KM,0:KM), . F(3,0:KM,0:KM), D1F(3,0:KM,0:KM), . D2F(3,0:KM,0:KM), SA1(KM,KM), . SA2(KM,KM), SAMIN COMMON/PHCOM1/ D, K COMMON/PHCOM2/ F0 COMMON/PHCOM3/ F COMMON/PHCOM4/ D1F COMMON/PHCOM5/ D2F COMMON/PHCOM6/ SA1, SA2, SAMIN C C Compute F = F0 + H*D. C DO 2 J = 0,K DO 1 I = 0,K F(1,I,J) = F0(1,I,J) + H*D(1,I,J) F(2,I,J) = F0(2,I,J) + H*D(2,I,J) F(3,I,J) = F0(3,I,J) + H*D(3,I,J) 1 CONTINUE 2 CONTINUE C C Compute PHI(F). C PHVAL = PHI (KM,K,F, D1F,D2F,SA1,SA2,SAMIN) RETURN END SUBROUTINE RES1 (A,B,C,D,E,S, R ) DOUBLE PRECISION A(3), B(3), C(3), D(3), E(3), S, R(3) 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 05/01/93 C C Given vectors A, B, C, D, E, and R of length 3 and a C scalar S, this subroutine overwrites R with C C R + S*(A X [(B X C) - (D X E)]) , C C where A X B denotes the vector cross product. The oper- C ation count is 15 adds and 21 multiplies. C C Modules required by RES1: None C C*********************************************************** C DOUBLE PRECISION F(3), G(3) C C F = B X C C F(1) = B(2)*C(3) - B(3)*C(2) F(2) = B(3)*C(1) - B(1)*C(3) F(3) = B(1)*C(2) - B(2)*C(1) C C G = F - (D X E) C G(1) = F(1) - D(2)*E(3) + D(3)*E(2) G(2) = F(2) - D(3)*E(1) + D(1)*E(3) G(3) = F(3) - D(1)*E(2) + D(2)*E(1) C C R = R + S*(A X G) C R(1) = R(1) + S*(A(2)*G(3)-A(3)*G(2)) R(2) = R(2) + S*(A(3)*G(1)-A(1)*G(3)) R(3) = R(3) + S*(A(1)*G(2)-A(2)*G(1)) RETURN END SUBROUTINE RES2 (A,B,C,D,E,F,S, R,T ) DOUBLE PRECISION A(3), B(3), C(3), D(3), E(3), F(3), . S, R(3), T(6) 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 03/29/94 C C This routine updates the order-3 matrix T and residual R C associated with a block of the linear system defined in C Subroutine SOR2. Given vectors A, B, C, D, E, F, and R of C length 3, an order-3 matrix T, and a scalar S, R and T are C overwritten with C C R + S*[A X G + H X F + *P] and C C T + S*(*I - A*A**T - P*P**T) , C C respectively, where A X B denotes vector cross product, C = (A**T)*A, A**T denotes the transpose of A, I C denotes the identity matrix, C C G = (B X C) - (D X E) , C H = C X E , and C P = S*H X A . C C The operation count is 47 adds and 66 multiplies. C C For S = 1/Abs(C X E), abs(S*H) = 1, and the matrix added C to T is symmetric and positive semi-definite. A compact C storage scheme is used with T stored as a vector of length C 6 containing the lower triangle stored by columns, or C equivalently, the upper triangle stored by rows: C (T11, T21, T31, T22, T32, T33). C C Modules required by RES2: None C C*********************************************************** C DOUBLE PRECISION A1S, A2S, A3S, G(3), H(3), P(3), . Q(3), QTG C C G = (B X C) - (D X E) C G(1) = B(2)*C(3) - B(3)*C(2) - D(2)*E(3) + D(3)*E(2) G(2) = B(3)*C(1) - B(1)*C(3) - D(3)*E(1) + D(1)*E(3) G(3) = B(1)*C(2) - B(2)*C(1) - D(1)*E(2) + D(2)*E(1) C C H = C X E C H(1) = C(2)*E(3) - C(3)*E(2) H(2) = C(3)*E(1) - C(1)*E(3) H(3) = C(1)*E(2) - C(2)*E(1) C C Q = S*H, P = Q X A, and QTG = 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