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 functional
C
C       PHI(F) = I( Abs(D1(F) X D2(F))**2 )/2 ,
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.  Function PHI also
C returns the surface area
C
C       SA = I( Abs(D1(F) X D2(F)) ) .
C
C Assuming the minimal surface with the given boundary con-
C ditions is unique, the minimizer of SA is unique only up
C to the parameterization, but the minimizer of PHI is the
C unique minimal surface (minimizer of SA) associated with a
C uniform parameterization -- constant Abs(D1(F) X D2(F)).
C The L2-gradient of SA is proportional to the negative
C mean curvature vector, and the minimal surface is thus
C characterized by zero mean curvature.  In place of the
C standard L2-gradient of PHI, the method employs an
C F-gradient G(F) corresponding to an inner product that de-
C pends on the evaluation point F.  It is thus a variable
C metric method.  Refer to Subroutine SOR1 for a description
C of the Sobolev gradient and SOR2 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 The nearly uniform parameterization, along with the uni-
C form grid, results in nearly equal-area triangular surface
C facets.
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 168*(KM+1)**2 = 437,000 for
C      KM = 50.  Note that a change in KM must also be made
C      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 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 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 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, SA, 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, SA
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 must not be smaller than the machine
C   precision.
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,SA)
C
C Compute the negative L2-gradient 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, 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) SA, PHIF, NIT, DMAX, RMSC, RMSG
      IF (INTER) THEN
        WRITE (*,180) NG
        WRITE (*,184) SA, 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:               SA = ',
     .           D12.6/
     .        5X,'Functional:                  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 L2-gradient:  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 H, SA, PHIF, and DPHI.
C
        WRITE (LPRT,190) H, SA, PHIF, DPHI
        IF (INTER) WRITE (*,190) H, SA, PHIF, DPHI
  190   FORMAT (/5X,'Step-size:',24X,'H = ',D9.3/
     .          5X,'Surface area of F+H*D:           SA = ',
     .             D12.6/
     .          5X,'Functional:              PHI(F+H*D) = ',
     .             D21.15/
     .          5X,'(PHI(F)-PHI(F+H*D))/PHI(F):    DPHI = ',
     .             D9.2/
     .          1X,'______________________________________',
     .             '___________________________')
      ENDIF
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 the negative L2-gradient 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, 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, 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) SA, PHIF, NIT, DMAX, RMSC, RMSG
        IF (INTER) WRITE (*,184) SA, 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 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) SA, PHI0, NIT, DMAX, RMSC, RMSG
        IF (INTER) WRITE (*,184) SA, PHI0, NIT, DMAX, RMSC,
     .                           RMSG
        WRITE (LPRT,190) H, SA, PHIF, DPHI
        IF (INTER) WRITE (*,190) H, SA, 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,FMTOL,SEARCH, H, PHIF,IER)
      INTEGER IER
      LOGICAL SEARCH
      DOUBLE PRECISION PHI0, FMTOL, 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/26/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       FMTOL = 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,SA = (PHCOM parameters) Output values from a
C                    call to PHI 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, FMTOL < 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
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), SA
      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, SA
      EXTERNAL PHVAL
C
C Test for invalid input.
C
      IF (K .LT. 2  .OR.  K .GT. KM  .OR.  FMTOL .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,FMTOL)
      ENDIF
C
C Compute the final value PHI(F0+H*D).
C
      PHIF = PHVAL(H)
      IER = 0
      RETURN
      END
      SUBROUTINE MAT (A, C )
      DOUBLE PRECISION A(3), 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                                                   04/26/92
C
C   Given a vector A of length 3 and a symmetric matrix C,
C this subroutine overwrites C with
C
C       C + (<A,A>I - A*A**T) ,
C
C where <A,A> = (A**T)*A and A**T denotes the transpose of
C A.  The operation count is 9 adds and 6 multiples.
C
C   The matrix added to C is symmetric and positive semi-
C definite and proportional to the orthogonal projection
C onto the complement of A.  A compact storage scheme is
C used with C stored as a vector of length 6 containing the
C lower triangle stored by columns, or equivalently, the
C 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) + A2S + A3S
      C(2) = C(2) - A(2)*A(1)
      C(3) = C(3) - A(3)*A(1)
      C(4) = C(4) + A1S + A3S
      C(5) = C(5) - A(3)*A(2)
      C(6) = C(6) + A1S + A2S
      RETURN
      END
      DOUBLE PRECISION FUNCTION PHI (KM,K,F, D1F,D2F,SA)
      INTEGER KM, K
      DOUBLE PRECISION F(3,0:KM,0:KM), D1F(3,0:KM,0:KM),
     .                 D2F(3,0:KM,0:KM), SA
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/29/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, the area
C SA of the piecewise linear interpolant on the triangula-
C tion defined by partitioning each cell by the diagonal
C with slope -1, and the functional whose minimum also min-
C imizes SA and results in a uniform parameterization:
C
C       PHI(F) = (K**2/4)*Sum( Abs(D1(F) X D2(F))**2 ) ,
C
C where the sum is over the 2*K**2 triangles, and Abs(D1(F)
C X D2(F))/2 is the surface area associated with a triangle.
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       SA = Surface area of the piecewise linear interpo-
C            lant.
C
C       PHI = Functional value at F defined above, or -1
C             if K < 2 or K > KM, in which case D1F, D2F,
C             and SA are not defined.
C
C Modules required by PHI:  None
C
C Intrinsic functions called by PHI:  DBLE, DSQRT
C
C***********************************************************
C
      DOUBLE PRECISION FN1, FN2, FN3, SASQ, SUM, SUMSQ
      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 SUM (sum of surface triangle areas) and SUMSQ
C   (sum of squared surface triangle areas).
C
      SUM = 0.D0
      SUMSQ = 0.D0
C
C Store D1F(i,j) and D2F(i,j) for i,j = 1 to K, and
C   accumulate SUM and SUMSQ.
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)
          SASQ = FN1*FN1 + FN2*FN2 + FN3*FN3
          SUMSQ = SUMSQ + SASQ
          SUM = SUM + DSQRT(SASQ)
          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)
          SASQ = FN1*FN1 + FN2*FN2 + FN3*FN3
          SUMSQ = SUMSQ + SASQ
          SUM = SUM + DSQRT(SASQ)
    2     CONTINUE
    3   CONTINUE
C
C Store SA and PHI.
C
      SA = SUM/2.D0
      PHI = DBLE(KK*KK)*SUMSQ/4.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, and SA.
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, and SA are updated.
C
C Module required by PHVAL:  PHI
C
C Common blocks required by PHVAL:  PHCOM1, PHCOM2, PHCOM3,
C                                   PHCOM4, PHCOM5
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), SA
      COMMON/PHCOM1/ D, K
      COMMON/PHCOM2/ F0
      COMMON/PHCOM3/ F
      COMMON/PHCOM4/ D1F
      COMMON/PHCOM5/ D2F, SA
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,SA)
      RETURN
      END
      SUBROUTINE RES1 (A,B,C,D,E, R )
      DOUBLE PRECISION A(3), B(3), C(3), D(3), E(3), 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                                                   08/17/93
C
C   Given vectors A, B, C, D, E, and R of length 3, this
C subroutine overwrites R with
C
C       R + 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 18 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 + (A X G)
C
      R(1) = R(1) + A(2)*G(3)-A(3)*G(2)
      R(2) = R(2) + A(3)*G(1)-A(1)*G(3)
      R(3) = R(3) + A(1)*G(2)-A(2)*G(1)
      RETURN
      END
      SUBROUTINE RES2 (A,B,C, R )
      DOUBLE PRECISION A(3), B(3), C(3), 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                                                   03/25/94
C
C   Given vectors A, B, C, and R of length 3, this subrou-
C tine overwrites R with
C
C       R + (A X B) X C ,
C
C where A X B denotes the vector cross product.  The oper-
C ation count is 9 adds and 12 multiplies.
C
C Modules required by RES2:  None
C
C***********************************************************
C
      DOUBLE PRECISION D(3)
C
C D = A X B
C
      D(1) = A(2)*B(3) - A(3)*B(2)
      D(2) = A(3)*B(1) - A(1)*B(3)
      D(3) = A(1)*B(2) - A(2)*B(1)
C
C R = R + (D X C)
C
      R(1) = R(1) + D(2)*C(3)-D(3)*C(2)
      R(2) = R(2) + D(3)*C(1)-D(1)*C(3)
      R(3) = R(3) + D(1)*C(2)-D(2)*C(1)
      RETURN
      END
      SUBROUTINE SOR1 (KM,K,F,D1F,D2F,OMEGA, NIT,DMAX, 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,
     .                 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/02/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) ] +
C          .5*D2[ D1(F) X D2(G) X D1(F) +
C                (D1(G) X D2(F)) X D1(F) ] = L(F)  and
C
C   G = 0 on the boundary of the unit square,
C
C where D1 and D2 are first partial derivative operators and
C L(F) is the negative L2-gradient of PHI(F).  All vectors,
C F, G, D1(F), etc. have length 3.
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       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.
C
C       C = Array containing grid-point values of the nega-
C           tive L2-gradient with zeros on the boundary.
C           C is not altered if IER = -1.
C
C       CTC = Sum (over the grid points) of squares of the
C             Euclidean norms of the grid-point L2-gradient
C             values unless IER = -1.
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.
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
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 L2-gradients 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 Copy F into G as the initial estimate of U = F-G, and
C   store zero boundary values in C.
C
      DO 2 J = 0,KK
        DO 1 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
    1     CONTINUE
    2   CONTINUE
C
C Top of iteration loop:  initialize DMX, and loop on
C                         interior grid points.
C
    3 DMX = 0.D0
      DO 12 J = 1,KM1
        JM1 = J-1
        JP1 = J+1
        DO 11 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 4 L = 1,3
            A(L) = 0.D0
            A(L+3) = 0.D0
            R(L) = 0.D0
            U(L) = G(L,I,J)
    4       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
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
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
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
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
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
C
          DO 5 L = 1,3
            B(L) = G(L,IP1,J) - U(L)
            D(L) = G(L,IP1,J) - G(L,IP1,JM1)
    5       CONTINUE
          CALL MAT (D2F(1,IP1,J), A )
          CALL RES1 (D2F(1,IP1,J),B,D2F(1,IP1,J),D,
     .               D1F(1,IP1,J), R )
C
          DO 6 L = 1,3
            T(L) = D2F(L,I,JP1) - D1F(L,IP1,J)
            D(L) = G(L,I,JP1) - U(L)
    6       CONTINUE
          CALL MAT (T, A )
          CALL RES1 (T,B,D2F(1,I,JP1),D,
     .               D1F(1,IP1,J), R )
C
          DO 7 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)
    7       CONTINUE
          CALL MAT (T, A )
          CALL RES1 (T,B,D2F(1,I,J),D,
     .               D1F(1,I,J), R )
C
          DO 8 L = 1,3
            D(L) = G(L,IM1,J) - G(L,IM1,JP1)
    8       CONTINUE
          CALL MAT (D2F(1,IM1,JP1), A )
          CALL RES1 (D2F(1,IM1,JP1),B,D2F(1,IM1,JP1),D,
     .               D1F(1,I,J), R )
C
          DO 9 L = 1,3
            B(L) = G(L,I,JP1) - U(L)
            D(L) = G(L,I,JP1) - G(L,IM1,JP1)
    9       CONTINUE
          CALL MAT (D1F(1,I,JP1), A )
          CALL RES1 (D1F(1,I,JP1),B,D1F(1,I,JP1),D,
     .               D2F(1,I,JP1), R )
C
          DO 10 L = 1,3
            B(L) = G(L,I,JM1) - U(L)
            D(L) = G(L,I,JM1) - G(L,IP1,JM1)
   10       CONTINUE
          CALL MAT (D1F(1,IP1,JM1), A )
          CALL RES1 (D1F(1,IP1,JM1),B,D1F(1,IP1,JM1),D,
     .               D2F(1,I,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 11
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)
   11     CONTINUE
   12   CONTINUE
      IF (ITMAX .EQ. 0) THEN
C
C Return G = -C for ITMAX = 0.
C
        DO 14 J = 0,KK
          DO 13 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)
   13       CONTINUE
   14     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 3
C
C Compute G = F-U and GSUM = GTG.
C
        GSUM = 0.D0
        DO 16 J = 0,KK
          DO 15 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
   15       CONTINUE
            IF (GSUM .GT. GTGMAX) GO TO 33
   16     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 -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
      END
      SUBROUTINE SOR2 (KM,K,F,D1F,D2F,OMEGA, NIT,DMAX, 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,
     .                 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) ]/3 +
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) ]/3 = L(F)  and
C
C   G = 0 on the boundary of the unit square,
C
C where D1 and D2 are first partial derivative operators and
C L(F) is the negative L2-gradient of PHI(F).  All vectors,
C F, G, D1(F), etc. have length 3.
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       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.
C
C       C = Array containing grid-point values of the nega-
C           tive L2-gradient with zeros on the boundary.
C           C is not altered if IER = -1.
C
C       CTC = Sum (over the grid points) of squares of the
C             Euclidean norms of the grid-point L2-gradient
C             values unless IER = -1.
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.
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
C MINSURF modules required by SOR2:  MAT, RES1, 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 L2-gradients 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 Copy F into G as the initial estimate of U = F-G,
C   store zero boundary values in C, and initialize
C   the residual norm squared RTR.
C
      DO 2 J = 0,KK
        DO 1 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
    1     CONTINUE
    2   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
    3 RTR0 = RTR
      RTR = 0.D0
      DMX = 0.D0
      DO 12 J = 1,KM1
        JM1 = J-1
        JP1 = J+1
        DO 11 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 4 L = 1,3
            A(L) = 0.D0
            A(L+3) = 0.D0
            R(L) = 0.D0
            U(L) = G(L,I,J)
    4       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    +[D2F(i+1,j) X D1F(i+1,j)] X (U(i+1,j)-U(i+1,j-1)) }
C
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    +[D2F(i,j+1) X D1F(i+1,j)] X (U(i,j+1)-U(i+1,j)) }
C
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    +[D2F(i,j) X D1F(i,j)] X (U(i-1,j)-U(i,j-1)) }
C
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    +[D2F(i-1,j+1) X D1F(i,j)] X (U(i-1,j+1)-U(i-1,j)) }
C
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    +[D1F(i,j+1) X D2F(i,j+1)] X (U(i,j+1)-U(i-1,j+1)) }
C
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    +[D1F(i+1,j-1) X D2F(i,j)] X (U(i+1,j-1)-U(i,j-1)) }
C
C
          DO 5 L = 1,3
            B(L) = G(L,IP1,J) - U(L)
            D(L) = G(L,IP1,J) - G(L,IP1,JM1)
    5       CONTINUE
          CALL MAT (D2F(1,IP1,J), A )
          CALL RES1 (D2F(1,IP1,J),B,D2F(1,IP1,J),D,
     .               D1F(1,IP1,J), R )
          CALL RES2 (D2F(1,IP1,J),D1F(1,IP1,J),D, R )
C
          DO 6 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)
    6       CONTINUE
          CALL MAT (T, A )
          CALL RES1 (T,B,D2F(1,I,JP1),D,
     .               D1F(1,IP1,J), R )
          CALL RES2 (D2F(1,I,JP1),D1F(1,IP1,J),E, R )
C
          DO 7 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)
    7       CONTINUE
          CALL MAT (T, A )
          CALL RES1 (T,B,D2F(1,I,J),D,
     .               D1F(1,I,J), R )
          CALL RES2 (D2F(1,I,J),D1F(1,I,J),E, R )
C
          DO 8 L = 1,3
            D(L) = G(L,IM1,J) - G(L,IM1,JP1)
    8       CONTINUE
          CALL MAT (D2F(1,IM1,JP1), A )
          CALL RES1 (D2F(1,IM1,JP1),B,D2F(1,IM1,JP1),D,
     .               D1F(1,I,J), R )
          CALL RES2 (D2F(1,IM1,JP1),D1F(1,I,J),D, R )
C
          DO 9 L = 1,3
            B(L) = G(L,I,JP1) - U(L)
            D(L) = G(L,I,JP1) - G(L,IM1,JP1)
    9       CONTINUE
          CALL MAT (D1F(1,I,JP1), A )
          CALL RES1 (D1F(1,I,JP1),B,D1F(1,I,JP1),D,
     .               D2F(1,I,JP1), R )
          CALL RES2 (D1F(1,I,JP1),D2F(1,I,JP1),D, R )
C
          DO 10 L = 1,3
            B(L) = G(L,I,JM1) - U(L)
            D(L) = G(L,I,JM1) - G(L,IP1,JM1)
   10       CONTINUE
          CALL MAT (D1F(1,IP1,JM1), A )
          CALL RES1 (D1F(1,IP1,JM1),B,D1F(1,IP1,JM1),D,
     .               D2F(1,I,J), R )
          CALL RES2 (D1F(1,IP1,JM1),D2F(1,I,J),D, 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)/3.D0
            C(2,I,J) = R(2)/3.D0
            C(3,I,J) = R(3)/3.D0
            CSUM = CSUM + R(1)**2 + R(2)**2 + R(3)**2
          ENDIF
          IF (ITMAX .EQ. 0) GO TO 11
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)
   11     CONTINUE
   12   CONTINUE
      IF (ITMAX .EQ. 0) THEN
C
C Return G = -C for ITMAX = 0.
C
        DO 14 J = 0,KK
          DO 13 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)
   13       CONTINUE
   14     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 17
        IF (DMX .GT. TOL  .AND.  ITER .LT. ITMAX) GO TO 3
C
C Compute G = (F-U)/3 and GSUM = GTG.
C
        GSUM = 0.D0
        DO 16 J = 0,KK
          DO 15 I = 0,KK
            G(1,I,J) = (F(1,I,J) - G(1,I,J))/3.D0
            G(2,I,J) = (F(2,I,J) - G(2,I,J))/3.D0
            G(3,I,J) = (F(3,I,J) - G(3,I,J))/3.D0
            GSUM = GSUM + G(1,I,J)**2 + G(2,I,J)**2 +
     .                    G(3,I,J)**2
   15       CONTINUE
            IF (GSUM .GT. GTGMAX) GO TO 33
   16     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
   17 NIT = ITER
      DMAX = DMX
      CTC = CSUM
      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
      END