dnsk.f

Go to the documentation of this file.
00001 C Work performed under the auspices of the U.S. Department of Energy
00002 C by Lawrence Livermore National Laboratory under contract number 
00003 C W-7405-Eng-48.
00004 C
00005       SUBROUTINE DNSK(X,Y,YPRIME,NEQ,RES,PSOL,WT,RPAR,IPAR,
00006      *   SAVR,DELTA,E,WM,IWM,CJ,SQRTN,RSQRTN,EPLIN,EPCON,
00007      *   S,CONFAC,TOLNEW,MULDEL,MAXIT,IRES,IERSL,IERNEW)
00008 C
00009 C***BEGIN PROLOGUE  DNSK
00010 C***REFER TO  DDASPK
00011 C***DATE WRITTEN   891219   (YYMMDD)
00012 C***REVISION DATE  900926   (YYMMDD)
00013 C***REVISION DATE  950126   (YYMMDD)
00014 C
00015 C
00016 C-----------------------------------------------------------------------
00017 C***DESCRIPTION
00018 C
00019 C     DNSK solves a nonlinear system of
00020 C     algebraic equations of the form
00021 C     G(X,Y,YPRIME) = 0 for the unknown Y.
00022 C
00023 C     The method used is a modified Newton scheme.
00024 C
00025 C     The parameters represent
00026 C
00027 C     X         -- Independent variable.
00028 C     Y         -- Solution vector.
00029 C     YPRIME    -- Derivative of solution vector.
00030 C     NEQ       -- Number of unknowns.
00031 C     RES       -- External user-supplied subroutine
00032 C                  to evaluate the residual.  See RES description
00033 C                  in DDASPK prologue.
00034 C     PSOL      -- External user-supplied routine to solve
00035 C                  a linear system using preconditioning. 
00036 C                  See explanation inside DDASPK.
00037 C     WT        -- Vector of weights for error criterion.
00038 C     RPAR,IPAR -- Real and integer arrays used for communication
00039 C                  between the calling program and external user
00040 C                  routines.  They are not altered within DASPK.
00041 C     SAVR      -- Work vector for DNSK of length NEQ.
00042 C     DELTA     -- Work vector for DNSK of length NEQ.
00043 C     E         -- Error accumulation vector for DNSK of length NEQ.
00044 C     WM,IWM    -- Real and integer arrays storing
00045 C                  matrix information such as the matrix
00046 C                  of partial derivatives, permutation
00047 C                  vector, and various other information.
00048 C     CJ        -- Parameter always proportional to 1/H (step size).
00049 C     SQRTN     -- Square root of NEQ.
00050 C     RSQRTN    -- reciprical of square root of NEQ.
00051 C     EPLIN     -- Tolerance for linear system solver.
00052 C     EPCON     -- Tolerance to test for convergence of the Newton
00053 C                  iteration.
00054 C     S         -- Used for error convergence tests.
00055 C                  In the Newton iteration: S = RATE/(1.D0-RATE),
00056 C                  where RATE is the estimated rate of convergence
00057 C                  of the Newton iteration.
00058 C
00059 C                  The closer RATE is to 0., the faster the Newton
00060 C                  iteration is converging; the closer RATE is to 1.,
00061 C                  the slower the Newton iteration is converging.
00062 C
00063 C                  The calling routine sends the initial value
00064 C                  of S to the Newton iteration.
00065 C     CONFAC    -- A residual scale factor to improve convergence.
00066 C     TOLNEW    -- Tolerance on the norm of Newton correction in
00067 C                  alternative Newton convergence test.
00068 C     MULDEL    -- A flag indicating whether or not to multiply
00069 C                  DELTA by CONFAC.
00070 C                  0  ==> do not scale DELTA by CONFAC.
00071 C                  1  ==> scale DELTA by CONFAC.
00072 C     MAXIT     -- Maximum allowed number of Newton iterations.
00073 C     IRES      -- Error flag returned from RES.  See RES description
00074 C                  in DDASPK prologue.  If IRES = -1, then IERNEW
00075 C                  will be set to 1.
00076 C                  If IRES < -1, then IERNEW will be set to -1.
00077 C     IERSL     -- Error flag for linear system solver.
00078 C                  See IERSL description in subroutine DSLVK.
00079 C                  If IERSL = 1, then IERNEW will be set to 1.
00080 C                  If IERSL < 0, then IERNEW will be set to -1.
00081 C     IERNEW    -- Error flag for Newton iteration.
00082 C                   0  ==> Newton iteration converged.
00083 C                   1  ==> recoverable error inside Newton iteration.
00084 C                  -1  ==> unrecoverable error inside Newton iteration.
00085 C-----------------------------------------------------------------------
00086 C
00087 C***ROUTINES CALLED
00088 C   RES, DSLVK, DDWNRM
00089 C
00090 C***END PROLOGUE  DNSK
00091 C
00092 C
00093       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
00094       DIMENSION Y(*),YPRIME(*),WT(*),DELTA(*),E(*),SAVR(*)
00095       DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*)
00096       EXTERNAL  RES, PSOL
00097 C
00098       PARAMETER (LNNI=19, LNRE=12)
00099 C
00100 C     Initialize Newton counter M and accumulation vector E.
00101 C
00102       M = 0
00103       DO 100 I=1,NEQ
00104 100     E(I) = 0.0D0
00105 C
00106 C     Corrector loop.
00107 C
00108 300   CONTINUE
00109       IWM(LNNI) = IWM(LNNI) + 1
00110 C
00111 C     If necessary, multiply residual by convergence factor.
00112 C
00113       IF (MULDEL .EQ. 1) THEN
00114         DO 320 I = 1,NEQ
00115 320       DELTA(I) = DELTA(I) * CONFAC
00116         ENDIF
00117 C
00118 C     Save residual in SAVR.
00119 C
00120       DO 340 I = 1,NEQ
00121 340     SAVR(I) = DELTA(I)
00122 C
00123 C     Compute a new iterate.  Store the correction in DELTA.
00124 C
00125       CALL DSLVK (NEQ, Y, X, YPRIME, SAVR, DELTA, WT, WM, IWM,
00126      *   RES, IRES, PSOL, IERSL, CJ, EPLIN, SQRTN, RSQRTN, RHOK,
00127      *   RPAR, IPAR)
00128       IF (IRES .NE. 0 .OR. IERSL .NE. 0) GO TO 380
00129 C
00130 C     Update Y, E, and YPRIME.
00131 C
00132       DO 360 I=1,NEQ
00133          Y(I) = Y(I) - DELTA(I)
00134          E(I) = E(I) - DELTA(I)
00135 360      YPRIME(I) = YPRIME(I) - CJ*DELTA(I)
00136 C
00137 C     Test for convergence of the iteration.
00138 C
00139       DELNRM = DDWNRM(NEQ,DELTA,WT,RPAR,IPAR)
00140       IF (DELNRM .LE. TOLNEW) GO TO 370
00141       IF (M .EQ. 0) THEN
00142         OLDNRM = DELNRM
00143       ELSE
00144         RATE = (DELNRM/OLDNRM)**(1.0D0/M)
00145         IF (RATE .GT. 0.9D0) GO TO 380
00146         S = RATE/(1.0D0 - RATE)
00147       ENDIF
00148       IF (S*DELNRM .LE. EPCON) GO TO 370
00149 C
00150 C     The corrector has not yet converged.  Update M and test whether
00151 C     the maximum number of iterations have been tried.
00152 C
00153       M = M + 1
00154       IF (M .GE. MAXIT) GO TO 380
00155 C
00156 C     Evaluate the residual, and go back to do another iteration.
00157 C
00158       IWM(LNRE) = IWM(LNRE) + 1
00159       CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR)
00160       IF (IRES .LT. 0) GO TO 380
00161       GO TO 300
00162 C
00163 C     The iteration has converged.
00164 C
00165 370    RETURN
00166 C
00167 C     The iteration has not converged.  Set IERNEW appropriately.
00168 C
00169 380   CONTINUE
00170       IF (IRES .LE. -2 .OR. IERSL .LT. 0) THEN
00171          IERNEW = -1
00172       ELSE
00173          IERNEW = 1
00174       ENDIF
00175       RETURN
00176 C
00177 C
00178 C------END OF SUBROUTINE DNSK-------------------------------------------
00179       END
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines