SUBROUTINE NLEQ1E(N,X,RTOL,IERR) C* Begin Prologue NLEQ1E INTEGER N DOUBLE PRECISION X(N), RTOL INTEGER IERR C ------------------------------------------------------------ C C* Title C C Numerical solution of nonlinear (NL) equations (EQ) C especially designed for numerically sensitive problems. C (E)asy-to-use driver routine for NLEQ1. C C* Written by U. Nowak, L. Weimann C* Purpose Solution of systems of highly nonlinear equations C* Method Damped affine invariant Newton method C (see references below) C* Category F2a. - Systems of nonlinear equations C* Keywords Nonlinear equations, Newton methods C* Version 2.3 C* Revision November 1991 C* Latest Change November 1991 C* Library CodeLib C* Code Fortran 77, Double Precision C* Environment Standard Fortran 77 environment on PC's, C workstations and hosts. C* Copyright (c) Konrad-Zuse-Zentrum fuer C Informationstechnik Berlin (ZIB) C Takustrasse 7, D-14195 Berlin-Dahlem C phone : + 49/30/84185-0 C fax : + 49/30/84185-125 C* Contact Bodo Erdmann C ZIB, Division Scientific Computing, C Department Numerical Analysis and Modelling C phone : + 49/30/84185-185 C fax : + 49/30/84185-107 C e-mail: erdmann@zib.de C C* References: C C /1/ P. Deuflhard: C Newton Methods for Nonlinear Problems. - C Affine Invariance and Adaptive Algorithms. C Series Computational Mathematics 35, Springer (2004) C C /2/ U. Nowak, L. Weimann: C A Family of Newton Codes for Systems of Highly Nonlinear C Equations - Algorithm, Implementation, Application. C ZIB, Technical Report TR 90-10 (December 1990) C C --------------------------------------------------------------- C C* Licence C You may use or modify this code for your own non commercial C purposes for an unlimited time. C In any case you should not deliver this code without a special C permission of ZIB. C In case you intend to use the code commercially, we oblige you C to sign an according licence agreement with ZIB. C C* Warranty C This code has been tested up to a certain level. Defects and C weaknesses, which may be included in the code, do not establish C any warranties by ZIB. ZIB does not take over any liabilities C which may follow from aquisition or application of this code. C C* Software status C This code is under care of ZIB and belongs to ZIB software class 1. C C ------------------------------------------------------------ C C* Summary: C ======== C Damped Newton-algorithm for systems of highly nonlinear C equations - damping strategy due to Ref. (1). C C (The iteration is done by subroutine N1INT actually. NLEQ1E C calls the standard interface driver NLEQ1, which itself does C some house keeping and builds up workspace.) C C Jacobian approximation by numerical differences. C C The numerical solution of the arising linear equations is C done by means of the subroutines *GEFA and *GESL ( Gauss- C algorithm with column-pivoting and row-interchange; C replace '*' by 'S' or 'D' for single or double precision C version respectively). C C ------------------------------------------------------------ C C* External subroutine to be supplied by the user C ============================================== C C FCN(N,X,F,IFAIL) Ext Function subroutine - the problem function C (The routine must be named exactly FCN !) C N Int Number of vector components (input) C Must not be altered! C X(N) Dble Vector of unknowns (input) C Must not be altered! C F(N) Dble Vector of function values (output) C IFAIL Int FCN evaluation-failure indicator. (output) C On input: Has always value 0 (zero). C On output: Indicates failure of FCN eval- C uation, if having a nonzero value. C If <0: NLEQ1E will be terminated with C IFAIL returned via IERR. C If =1: A new trial Newton iterate will be C computed, with the damping factor C reduced to it's half. C If =2: A new trial Newton iterate will C computed, with the damping factor C reduced by a reduction factor, C which must be output through F(1) C by FCN, and it's value must be C >0 and < 1. C C* Parameters list description C =========================== C C Name Type In/Out Meaning C C N Int In Number of unknowns ( N .LE. 50 ) C X(N) Dble In Initial estimate of the solution C Out Solution values ( or final values, C respectively ) C RTOL Dble In Required relative precision of C solution components - C RTOL.GE.EPMACH*TEN*N C Out Finally achieved accuracy C IERR Int Out Return value parameter C < 0 Termination forced by user function FCN C due to IFAIL<0 on output, IERR is set C to IFAIL C = 0 successfull completion of the iteration, C solution has been computed C > 0 see list of error/warning messages below C C* Error and warning messages: C =========================== C C 1 Termination, since Jacobian matrix became singular C 2 Termination after 100 iterations C 3 Termination, since damping factor became to small C 4 Warning: Superlinear or quadratic convergence slowed down C near the solution. C Iteration has been stopped therefore with an approximation C of the solution not such accurate as requested by RTOL, C because possibly the RTOL requirement may be too stringent C (i.e. the nonlinear problem is ill-conditioned) C 5 Warning: Iteration stopped with termination criterion C (using RTOL as requested precision) satisfied, but no C superlinear or quadratic convergence has been indicated yet. C Therefore, possibly the error estimate for the solution may C not match good enough the really achieved accuracy. C 20 Bad input value to parameter N; 1 .LE. N .LE. 50 required C 21 Nonpositive value for RTOL supplied C 82 Termination, because user routine FCN returned with IFAIL>0 C C Note : in case of failure: C - use better initial guess C - or refine model C - or use non-standard options and/or analytical Jacobian C via the standard interface NLEQ1 C C* Machine dependent constants used: C ================================= C C DOUBLE PRECISION EPMACH in N1PCHK, N1INT C DOUBLE PRECISION GREAT in N1PCHK C DOUBLE PRECISION SMALL in N1PCHK, N1INT, N1SCAL C C* Subroutines called: NLEQ1 C C ------------------------------------------------------------ C* End Prologue INTEGER NMAX, LIOPT, LIWK, LRWK, LUPRT PARAMETER ( NMAX=50 ) PARAMETER ( LIOPT=50, LIWK=NMAX+50, LRWK=(NMAX+13)*NMAX+60 ) PARAMETER ( LUPRT=6 ) EXTERNAL FCN, NLEQ1 DOUBLE PRECISION XSCAL(NMAX) INTEGER IOPT(LIOPT), IWK(LIWK) DOUBLE PRECISION RWK(LRWK) INTEGER I, NIW, NRW CHARACTER CHGDAT*20, PRODCT*8 C C Version: 2.3 Latest change: C ----------------------------------------- C DATA CHGDAT /'November 15, 1991 '/ DATA PRODCT /'NLEQ1E '/ C* Begin NIW = N+50 NRW = (N+13)*N+60 C Checking dimensional parameter N IF ( N.LT.1 .OR. N.GT.NMAX ) THEN WRITE(LUPRT,1001) NMAX, N 1001 FORMAT(/,' Error: Bad input to parameter N supplied', $ /,8X,'choose 1 .LE. N .LE. ',I3, $ ' , your input is: N = ',I5) IERR = 20 RETURN ENDIF DO 10 I=1,LIOPT IOPT(I) = 0 10 CONTINUE DO 20 I=1,NIW IWK(I) = 0 20 CONTINUE DO 30 I=1,NRW RWK(I) = 0.0D0 30 CONTINUE DO 40 I=1,N XSCAL(I) = 0.0D0 40 CONTINUE C Print errors, warnings, monitor and time monitor C to standard output IOPT(11) = 3 IOPT(12) = LUPRT IOPT(13) = 3 IOPT(14) = LUPRT IOPT(19) = 1 IOPT(20) = LUPRT C Maximum number of Newton iterations IWK(31) = 100 C CALL NLEQ1(N,FCN,DUMMY,X,XSCAL,RTOL,IOPT,IERR,LIWK,IWK,LRWK,RWK) IF (IERR.EQ.82 .AND. IWK(23).LT.0) IERR=IWK(23) C End of subroutine NLEQ1E RETURN END