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