function [x,info,wk] = nleq1(fcn,x,xscal,rtol,iopt,par,wk)
%     ------------------------------------------------------------
%
%*  Title
%
%     Numerical solution of nonlinear (NL) equations (EQ)
%     especially designed for numerically sensitive problems.
%
%*  Written by        U. Nowak, L. Weimann 
%*  Purpose           Solution of systems of highly nonlinear equations
%*  Method            Damped affine invariant Newton method
%                     (see references below)
%*  Category          F2a. - Systems of nonlinear equations
%*  Keywords          Nonlinear equations, Newton methods
%*  Version           2.4
%*  Revision          November 2006
%*  Latest Change     June 2017
%*  Library           CodeLib
%*  Code              MatLab 6.0
%*  Environment       Standard MatLab environment on PCs,
%                     workstations and hosts.
%*  Copyright     (c) Konrad-Zuse-Zentrum fuer
%                     Informationstechnik Berlin (ZIB)
%                     Takustrasse 7, D-14195 Berlin-Dahlem
%                     phone : + 49/30/84185-0
%                     fax   : + 49/30/84185-125
%*  Contact           Bodo Erdmann
%                     ZIB, Division Scientific Computing, 
%                          Department Numerical Analysis and Modelling
%                     phone : + 49/30/84185-185
%                     fax   : + 49/30/84185-107
%                     e-mail: erdmann@zib.de
%
%*    References:
%
%     /1/ P. Deuflhard:
%         Newton Methods for Nonlinear Problems. -
%         Affine Invariance and Adaptive Algorithms.
%         Series Computational Mathematics 35, Springer (2004)
%
%     /2/ U. Nowak, L. Weimann:
%         A Family of Newton Codes for Systems of Highly Nonlinear
%         Equations - Algorithm, Implementation, Application.
%         ZIB, Technical Report TR 90-10 (December 1990)
%
%  ---------------------------------------------------------------
%
%* Licence
%    You may use or modify this code for your own non commercial
%    purposes for an unlimited time. 
%    In any case you should not deliver this code without a special 
%    permission of ZIB.
%    In case you intend to use the code commercially, we oblige you
%    to sign an according licence agreement with ZIB.
%
%* Warranty 
%    This code has been tested up to a certain level. Defects and
%    weaknesses, which may be included in the code, do not establish
%    any warranties by ZIB. ZIB does not take over any liabilities
%    which may follow from acquisition or application of this code.
%
%* Software status 
%    This code is under care of ZIB and belongs to ZIB software class 1.
%
%     ------------------------------------------------------------
%
%*    Summary:
%     ========
%     Damped Newton-algorithm for systems of highly nonlinear
%     equations - damping strategy due to Ref. (1).
%
%     (The iteration is done by function N1INT currently. NLEQ1
%      itself does some house keeping and builds up workspace.)
%
%     Jacobian approximation by numerical differences or user
%     supplied function FCN.
%
%     The numerical solution of the arising linear equations is
%     done by means of MatLab builtin lu-factorization routine in the
%     dense or sparse matrix case; or by the functions DGBFA and DGBSL,
%     which have been converted from the equal named Fortran LINPACK 
%     subroutines into MatLab code, in the band matrix case.
%     For special purposes these routines may be substituted.
%
%     This is a driver routine for the core solver N1INT.
%
%     ------------------------------------------------------------
%
%*    Parameters list description (* marks inout parameters)
%     ======================================================
%
%*    External functions (to be supplied by the user)
%     ===============================================
% 
%     [FOUT,FAIL] =FCN(X,FLAG,PAR) String  Name of problem- and Jacobian function
%       X()           Float   Vector of unknowns (input)
%       FLAG          String  Operation flag - the following values must be
%                             supported:
%                             '' (empty string) : FOUT must return the problem-
%                                                 function value f(X) ;
%                             'jacobian' : FOUT must return the associated 
%                                          Jacobian Jac(x)
%       PAR            AnyType  A (required!) user parameter
%
%       FOUT()         Float  Vector of returned function values or Jacobian matrix
%       FAIL           Int    evaluation-failure indicator. (output)
%                             On output: Indicates failure of FCN eval-
%                                uation, if having a value <= 2.
%                             If <0: NLEQ1 will be terminated with 
%                                    error code = 82, and FAIL stored
%                                    to wk.ifail.
%                             If =1: A new trial Newton iterate will
%                                    computed, with the damping factor
%                                    reduced to its half.
%                             If =2: A new trial Newton iterate will
%                                    computed, with the damping factor
%                                    reduced by a reduct. factor, which
%                                    must be output through F(1) by FCN,
%                                    and it is value must be >0 and < 1.
%                             Note, that if FAIL = 1 or 2, additional
%                             conditions concerning the damping factor,
%                             e.g. the minimum damping factor or the
%                             bounded damping strategy may also influ-
%                             ence the value of the reduced damping 
%                             factor.
%
%*    Input parameters of NLEQ1
%     =========================
%
%     FCN            String Name of problem- and Jacobian function (see above)
%     X(1:N)         Dble   Initial estimate of the solution
%     XSCAL(1:N)     Dble   User scaling (lower threshold) of the 
%                           iteration vector X(N)
%     RTOL           Dble   Required relative precision of
%                           solution components -
%                           RTOL >= EPMACH*TEN*N
%     IOPT           Struct Structure of run-time options. Omit options
%                           to get the default values (details see below)
%     PAR            Any    A user function parameter, which is passed on
%                           every call of the user function as the third
%                           input argument.
%     WK             Struct A structured workspace. Normally an empty
%                           array should be passed to this argument.
%                           However, on successor steps in onestep mode,
%                           the output parameter WK from the previous step
%                           call must be passed to this parameter.
%                           
%
%*    Output parameters of NLEQ1
%     ==========================
%
%     X(N)            Dble   Solution values ( or final values,
%                            respectively )
%     INFO            Struct A structure of additional output values.
%                            The structure has the following fields:
%        XSCAL(N)     Dble   After return with IERR >= 0, it contains
%                            the latest internal scaling vector used
%                            After return with IERR == -1 in onestep-
%                            mode it contains a possibly adapted 
%                            (as described below) user scaling vector:
%                            If (XSCAL(I) <  SMALL) XSCAL(I) = SMALL ,
%                            If (XSCAL(I) >  GREAT) XSCAL(I) = GREAT .
%                            For SMALL and GREAT, see section machine
%                            constants below  and regard note 1.
%        RTOL         Dble   Finally achieved (relative) accuracy.
%                            The estimated absolute error of component i
%                            of x_out is approximately given by
%                             abs_err(i) = RTOL * XSCAL_out(i) ,
%                            where (approximately)
%                             XSCAL_out(i) = 
%                                max(abs(X_out(i)),XSCAL_in(i)).
%        IERR         Int    Return value parameter
%                            =-1 sucessfull completion of one iteration
%                                step, subsequent iterations are needed 
%                                to get a solution. (stepwise mode only) 
%                            = 0 successfull completion of the iteration,
%                                solution has been computed
%                            > 0 see list of error messages below
%        XITER(N,:)   Dble   An array holding all iterates of the
%                            Newton iteration run, stored as column
%                            vectors. 
%        NATLEVEL(:)  Dble   The sequence of natural levels of the 
%                            newton corrections over the iteration steps
%        SIMLEVEL(:)  Dble   The sequence of natural levels of the simplified 
%                            newton corrections over the iteration steps 
%        STDLEVEL(:)  Dble   The sequence of standard levels over the
%                            iteration steps
%        PRECISION(:) Dble   The sequence of acheived precisions over
%                            the iteration steps.
%        DAMPINGFC(:) Dble   The sequence of accepted damping factors
%                            over the iteration steps.
%        NITER        Int    Number of Newton-iterations done
%        NCORR        Int    Number of corrector steps done
%        NREJR1       Int    Number of rejected Newton iteration steps
%                             done with a rank-1 approximated Jacobian
%        NJAC         Int    Number of Jacobian generations or JAC-calls done
%        NFCN         Int    Number of FCN-evaluations done
%        NFCNJ        Int    Number of FCN-evaluations for Jacobian
%                             approximation done
%
%     Note 1.
%        The machine dependent values SMALL, GREAT and EPMACH are
%        gained from calls of the machine constants function D1MACH.
%        As delivered, this function is adapted to use constants 
%        suitable for all machines with IEEE arithmetic. If you use
%        another type of machine, you have to change the DATA state-
%        ments for IEEE arithmetic in D1MACH into comments and to 
%        uncomment the set of DATA statements suitable for your machine.
%
%*   Options IOPT.Name
%    =================
%
%          Name   Default  Meaning
%
%          mode   0        =0 Standard mode initial call:
%                             Return when the required accuracy for the
%                             iteration vector is reached. User defined
%                             parameters are evaluated and checked.
%                             Standard mode successive call:
%                             If NLEQ1 was called previously with MODE=1,
%                             it performs all remaining iteration steps.
%                          =1 Stepwise mode:
%                             Return after one Newton iteration step.
%          jacgen 0        Method of Jacobian generation
%                          =0 Standard method is JACGEN=2
%                          =1 User supplied function FCN will be 
%                             called to generate Jacobian matrix
%                          =2 Jacobian approximation by numerical
%                             differentation (see functions N1JAC
%                             and N1JACB)
%                          =3 Jacobian approximation by numerical
%                             differentation with feedback control
%                             (see functions N1JCF and N1JCFB)
%          mstor  0        =0 The Jacobian A is a dense matrix or
%                             sparse matrix (note that in case of a
%                             sparse matrix, jacgen must be set to 1,
%                             and the sparse Jacobian must be supplied
%                             analytically by function FCN)
%                          =1 A is a band matrix
%          ml     0        Lower bandwidth of A (excluding the
%                          diagonal);
%                          IOPT.ml ignored, if IOPT.mstor ~= 1
%          mu     0        Upper bandwidth of A (excluding the
%                          diagonal);
%                          IOPT.mu ignored, if IOPT.mstor ~= 1
%          iscal  0        Determines how to scale the iterate-vector:
%                          =0 The user supplied scaling vector XSCAL is
%                             used as a (componentwise) lower threshold
%                             of the current scaling vector
%                          =1 The vector XSCAL is always used as the
%                             current scaling vector
%          mprerr 0        Print error messages
%                          =0 No output
%                          =1 Error messages
%                          =2 Warnings additionally
%                          =3 Informal messages additionally
%          luerr  1        Logical unit number for error messages
%          mprmon 0        Print iteration Monitor
%                          =0 No output
%                          =1 Standard output
%                          =2 Summary iteration monitor additionally
%                          =3 Detailed iteration monitor additionally
%                          =4,5,6 Outputs with increasing level addi-
%                             tional increasing information for code
%                             testing purposes. Level 6 produces
%                             in general extremely large output!
%          lumon  1        Logical unit number for iteration monitor
%          mprsol 0        Print solutions
%                          =0 No output
%                          =1 Initial values and solution values
%                          =2 Intermediate iterates additionally
%          lusol  1        Logical unit number for solutions
%          mprtim 0        Output level for the time monitor
%                          = 0 : no time measurement and no output
%                          = 1 : time measurement will be done and
%                                summary output will be written -
%                                regard note 5a.
%          lutim  1        Logical output unit for time monitor
%          nonlin 3        Problem type specification
%                          =1 Linear problem
%                             Warning: If specified, no check will be
%                             done, if the problem is really linear, and
%                             NLEQ1 terminates unconditionally after one
%                             Newton-iteration step.
%                          =2 Mildly nonlinear problem
%                          =3 Highly nonlinear problem
%                          =4 Extremely nonlinear problem
%          qrank1 0        =0 (.FALSE.) Rank-1 updates by Broyden-
%                             approximation are inhibited.
%                          =1 (.TRUE.) Rank-1 updates by Broyden-
%                             approximation are allowed.
%          qordi  0        =0 (.FALSE.) Standard program mode 
%                          =1 (.TRUE.)  Special program mode:
%                             Ordinary Newton iteration is done, e.g.:
%                             No damping strategy and no monotonicity
%                             test is applied
%          qsimpl 0        =0 (.FALSE.) Standard program mode
%                          =1 (.TRUE.)  Special program mode:
%                             Simplified Newton iteration is done, e.g.:
%                             The Jacobian computed at the starting
%                             point is fixed for all subsequent
%                             iteration steps, and
%                             no damping strategy and no monotonicity
%                             test is applied.
%          norowscal 0     Inhibit automatic row scaling: 
%                          =0 (.FALSE.) Automatic row scaling of
%                             the linear system is activ: 
%                             Rows i=1,...,N will be divided by
%                             max j=1,...,N (abs(a(i,j))) 
%                          =1 (.TRUE.) No row scaling of the linear
%                             system. Recommended only for well row-
%                             scaled nonlinear systems.
%          boundeddamp     Bounded damping strategy switch:
%                          =0 The default switch takes place, dependent
%                             on the setting of nonlin (=iopt.nonlin):
%                             nonlin = 0,1,2,3 -> boundeddamp = off ,
%                             nonlin = 4 -> boundeddamp = on
%                          =1 means always boundeddamp = on 
%                          =2 means always boundeddamp = off
%          iormon          Convergence order monitor 
%                          =0 Standard option is iormon=2 
%                          =1 Convergence order is not checked,
%                             the iteration will be always proceeded
%                             until the solution has the required 
%                             precision RTOL (or some error condition
%                             occured)
%                          =2 Use additional 'weak stop' criterion:
%                             Convergence order is monitored
%                             and termination due to slowdown of the
%                             convergence may occur.
%                          =3 Use additional 'hard stop' criterion:
%                             Convergence order is monitored
%                             and termination due to superlinear 
%                             convergence slowdown may occur. 
%                          In case of termination due to convergence
%                          slowdown, the warning code IERR=4 will be
%                          set.
%                          In cases, where the Newton iteration con-
%                          verges but superlinear convergence order has
%                          never been detected, the warning code IERR=5 
%                          is returned.
%
%     Note 3:
%         If NLEQ1 terminates with IERR=2 (maximum iterations)
%         or  IERR=3 (small damping factor), you may try to continue
%         the iteration by increasing wk.nitmax or decreasing wk.fcmin
%         (see WK) and setting wk.qsucc to 1.
%
%     Note 4 : Storage of user supplied banded Jacobian
%        In the band matrix case, the following lines may build
%        up the analytic Jacobian A;
%        Here AFL denotes the quadratic matrix A in dense form,
%        and ABD the rectangular matrix A in banded form :
%
%                   ml = iopt.ml ;
%                   mu = iopt.mu ;
%                   mh = mu+1
%                   for j = 1:n
%                     i1 = max(1,j-mu) ;
%                     i2 = min(n,j+ml) ;
%                     for i = i1:i2
%                       k = i-j+mh ;
%                       abd(k,j) = afl(i,j) ;
%                     end
%                   end
%
%         The total number of rows needed in  ABD  is  ML+MU+1 .
%         The  MU by MU  upper left triangle and the
%         ML by ML  lower right triangle are not referenced.
%
%
%*   Optional input/output in WK:
%    ============================
%
%          Name          Meaning
%
%          niter  IN/OUT Number of Newton-iterations
%          ncorr  IN/OUT Number of corrector steps
%          nfcn   IN/OUT Number of FCN-evaluations
%          njac   IN/OUT Number of Jacobian generations or JAC-calls
%          nfcnj  IN/OUT Number of FCN-evaluations for Jacobian
%                        approximation
%          nrejr1 IN/OUT Number of rejected Newton iteration steps
%                        done with a rank-1 approximated Jacobian
%          idcode IN/OUT Output: The 8 decimal digits program identi-
%                        fication number ppppvvvv, consisting of the
%                        program code pppp and the version code vvvv.
%                        Input: If containing a negative number,
%                        it will only be overwritten by the identi-
%                        fication number, immediately followed by
%                        a return to the calling program.      
%          ifail  OUT    Set in case of failure of N1FACT (IERR=80),
%                        N1SOLV (IERR=81), FCN (IERR=82) or JAC(IERR=83)
%                        to the nonzero IFAIL value returned by the 
%                        routine indicating the failure .
%          iconv  OUT    Current status of of the convergence monitor
%                        (only if convergence order monitor is on - 
%                         see IORMON(=IOPT(39)))
%                        =0: No convergence indicated yet
%                        =1: Damping factor is 1.0d0
%                        =2: Superlinear convergence in progress
%                        =3: Quadratic convergence in progress
%          nitmax IN     Maximum number of permitted iteration
%                        steps (default: 50)
%          new    IN/OUT Count of consecutive rank-1 updates
%          nbroy  IN     Maximum number of possible consecutive 
%                        iterative Broyden steps. Note that the total  
%                        size of workspace (wk.dxsave) needed depends 
%                        on this value.
%                        Default is N (see parameter N),
%                        if MSTOR=0 (=IOPT(4)), 
%                        and ML+MU+1 (=IOPT(6)+IOPT(7)+1), if MSTOR=1
%                        (but minimum is always 10) - 
%                        provided that Broyden steps are allowed. 
%                        If Broyden is inhibited, NBROY is always set to
%                        zero.
%          conv   OUT    The achieved relative accuracy after the  
%                        latest step
%          sumx   OUT    Natural level (not Normx of printouts)
%                        of the latest iterate, i.e. Sum(DX(i)**2),
%                        where DX = scaled Newton correction.
%          dlevf  OUT    Standard level (not Normf of printouts)
%                        of the latest iterate, i.e. Norm2(F(X)),
%                        where F =  nonlinear problem function.
%          fcbnd  IN     Bounded damping strategy restriction factor
%                        (Default is 10)
%          fcstrt IN     Damping factor for first Newton iteration -
%                        overrides option NONLIN, if set (see note 6)
%          fcmin  IN     Minimal allowed damping factor (see note 6)
%          sigma  IN     Broyden-approximation decision parameter
%                        Required choice: SIGMA >= 1. Increasing this
%                        parameter make it less probable that the algo-
%                        rith performs rank-1 updates.
%                        Rank-1 updates are inhibited, if 
%                        SIGMA > 1/FCMIN is set. (see note 6)
%          sigma2 IN     Decision parameter about increasing damping
%                        factor to corrector if predictor is small.
%                        Required choice: SIGMA2 >= 1. Increasing this
%                        parameter make it less probable that the algo-
%                        rith performs rank-1 updates.
%          ajdel  IN     Jacobian approximation without feedback:
%                        Relative pertubation for components
%                        (Default: sqrt(epmach*10), epmach: relative
%                         machine precision) 
%          ajmin  IN     Jacobian approximation without feedback:
%                        Threshold value (Default: 0.0d0)
%                          The absolute pertubation for component k is
%                          computed by 
%                          DELX := AJDEL*max(abs(Xk),AJMIN)
%         etadif  IN     Jacobian approximation with feedback:
%                        Target value for relative pertubation ETA of X
%                        (Default: 1.0d-6)
%         etaini  IN     Jacobian approximation with feedback:
%                        Initial value for denominator differences
%                        (Default: 1.0d-6)
%
%     Note 6:
%       The default values of the internal parameters may be obtained
%       from the monitor output with at least iopt.mprmon set to 2. 
%
%*   Error and warning messages:
%    ===========================
%
%      1    Termination, since jacobian matrix became singular
%      2    Termination after NITMAX iterations ( as indicated by
%           input parameter NITMAX=IWK(31) )
%      3    Termination, since damping factor became to small
%      4    Warning: Superlinear or quadratic convergence slowed down
%           near the solution.
%           Iteration has been stopped therefore with an approximation
%           of the solution not such accurate as requested by RTOL,
%           because possibly the RTOL requirement may be too stringent
%           (i.e. the nonlinear problem is ill-conditioned)
%      5    Warning: Iteration stopped with termination criterion 
%           (using RTOL as requested precision) satisfied, but no 
%           superlinear or quadratic convergence has been indicated yet.
%           Therefore, possibly the error estimate for the solution may
%           not match good enough the really achieved accuracy.
%     21    Nonpositive value for RTOL supplied
%     22    Negative scaling value via vector XSCAL supplied
%     30    One or more fields specified in IOPT are invalid
%           (for more information, see error-printout)
%     80    Error signalled by linear solver routine N1FACT,
%           for more detailed information see IFAIL-value
%           stored to IWK(23)
%     81    Error signalled by linear solver routine N1SOLV,
%           for more detailed information see IFAIL-value
%           stored to IWK(23)
%           (not used by standard routine N1SOLV)
%     82    Error signalled by user routine FCN (Nonzero value
%           returned via IFAIL-flag; stored to IWK(23) )
%     83    Error signalled by user routine JAC (Nonzero value
%           returned via IFAIL-flag; stored to IWK(23) )
%
%     Note 7 : in case of failure:
%        -    use non-standard options
%        -    or turn to Newton-algorithm with rank strategy
%        -    use another initial guess
%        -    or reformulate model
%        -    or apply continuation techniques
%
%*    Machine dependent constants used:
%     =================================
%
%     DOUBLE PRECISION EPMACH  in  N1PCHK, N1INT
%     DOUBLE PRECISION GREAT   in  N1PCHK
%     DOUBLE PRECISION SMALL   in  N1PCHK, N1INT, N1SCAL
%
%*    functions called: N1PCHK, N1INT
%
%     ------------------------------------------------------------
%*    End Prologue
%
%*    Summary of changes of the underlying Fortran code:
%     ==================================================
%      
%     2.2.1  91, May 30     Time monitor included
%     2.2.2  91, May 30     Bounded damping strategy implemented
%     2.2.3  91, June 19    AJDEL, AJMIN as RWK-options for JACGEN == 2,
%                           ETADIF, ETAINI as RWK-opts. for JACGEN == 3
%                           FCN-count changed for anal. Jacobian
%     2.2.4  91, August  9  Convergence order monitor included
%     2.2.5  91, August 13  Standard Broyden updates replaced by
%                           iterative Broyden
%     2.2.6  91, Sept.  16  Damping factor reduction by FCN-fail imple-
%                           mented
%     2.3    91, Dec.   20  New Release for CodeLib
%            92, March  11  Level of accepted simplified correction
%                           stored to RWK(IRWKI+4)
%            00, July   12  RTOL output-value bug fixed
%            06, Jan.   24  IERR=5 no longer returned if residuum of
%                           final iterate is exactly zero
%     2.4    06, Nov.   23  New code added to support best performance
%                           when user supplies a sparse Jacobian matrix,
%                           and other changes to generally enhance
%                           the MATLAB performance.
%            06, Nov.   24  Fixed missing info output parameter at early
%                           returns from NLEQ1 
%   
%     ------------------------------------------------------------
%
%     Further WK fields (only internally used)
%
%               Name     Meaning
%
%               FCKEEP   Damping factor of previous successfull iter.
%               FCA      Previous damping factor
%               FCPRI    A priori estimate of damping factor
%               DMYCOR   Number My of latest corrector damping factor
%                        (kept for use in rank-1 decision criterium)
%               SUMXS    natural level of accepted simplified correction
%
%     Internal arrays stored in WK (see routine N1INT for descriptions)
%
%               Array(size(s))  Type   Remarks
%
%               A(M1,N)         Perm   M1=N (full mode) or 
%                                    M1=2*iopt.ml+iopt.mu+1 (band mode)
%               DXSAVE(N,NBROY) Perm   NBROY=wk.nbroy (Default: N or 0)
%               DX(N)           Perm  
%               DXQ(N)          Perm 
%               XA(N)           Perm
%               F(N)            Perm
%               FW(N)           Perm
%               XWA(N)          Perm
%               FA(N)           Perm
%               ETA(N)          Perm   Only used for iopt.jacgen=3
%               XW(N)           Temp
%               DXQA(N)         Temp
%     
%
      persistent xiter sumxall dlevfall sumxqall tolall fcall ;
      iver=21112401 ;
%
%     Version: 2.4               Latest change: 
%     -----------------------------------------
%
      chgdat = 'November 24, 2006   ' ;
      prodct = 'nleq1   '             ;
%*    Begin
      ierr = 0 ;
      if isfield(wk,'iver')
        qvchk = wk.iver < 0 ;
        wk.iver = iver ;
        if qvchk
          info.ierr = 0 ;
          return ;
        end ;
      else
        wk.iver = iver ;
      end
      
%        Print error messages?
      [mprerr,iopt] = getopt(iopt,'mprerr',0) ;
      [luerr, iopt] = getopt(iopt,'luerr', 1) ;
      if luerr  ==  0
        luerr = 1 ;
        iopt.luerr = luerr ;
      end
%        print iteration monitor?
      [mprmon,iopt] = getopt(iopt,'mprmon',0) ;
      [lumon ,iopt] = getopt(iopt,'lumon', 1) ;
      if lumon  <=  0  ||  lumon  >  99
        lumon = 1 ;
        iopt.lumon = lumon ;
      end
%        print intermediate solutions?
      [mprsol,iopt] = getopt(iopt,'mprsol',0) ;
      [lusol, iopt] = getopt(iopt,'lusol', 1) ;
      if lusol  ==  0
        lusol = 1 ;
        iopt.lusol=lusol ;
      end
%        print time summary statistics?
      [mprtim,iopt] = getopt(iopt,'mprtim',0) ;
      [lutim, iopt] = getopt(iopt,'lutim', 1) ;
      if lutim  ==  0
        lutim = 1 ;
        iopt.lutim=lutim ;
      end
      [qsucc,wk] = getopt(wk,'qsucc',0) ;
      qinimo = mprmon >= 1 &&  ~ qsucc ;
%     Print NLEQ1 heading lines
      if qinimo
        fprintf(lumon,'%s\n\n%s\n\n','   N L E Q 1  *****  V e r s i o n  2 . 3 ***', ...
                      ' Newton-Method for the solution of nonlinear systems') ;
      end
%     Check input parameters and options
      n = length(x) ;
      [ierr,iopt,rtol,xscal] = n1pchk(n,x,xscal,rtol,iopt,wk) ;
%     Exit, if any parameter error was detected till here
      if ierr ~= 0
        info.ierr = ierr ;
        return ;
      end
%
      if ~ qsucc
        xiter    = [] ;
        sumxall  = [] ;
        dlevfall = [] ;
        sumxqall = [] ;
        tolall   = [] ;
        fcall    = [] ;
      end
%
      [mstor,iopt] = getopt(iopt,'mstor',0) ;
      if mstor == 0
        m1=n ;
        m2=n ;
      elseif mstor == 1
        [ml,iopt] = getopt(iopt,'ml',0) ;
        [mu,iopt] = getopt(iopt,'mu',0) ;
        m1=2*ml+mu+1 ;
        m2=ml+mu+1 ;
      end
      [jacgen,iopt] = getopt(iopt,'jacgen',0) ;
      if jacgen == 0
        jacgen=2 ;
      end
      iopt.jacgen = jacgen ;
      [qrank1,iopt] = getopt(iopt,'qrank1',0) ;
      [qordi,iopt]  = getopt(iopt,'qordi',0) ;
      [qsimpl,iopt] = getopt(iopt,'qsimpl',0) ;
      if qrank1
        [nbroy,wk] = getopt(wk,'nbroy',0) ;
        if nbroy == 0
          nbroy=max(m2,10) ;
        end
        wk.nbroy=nbroy ;
      else
        nbroy=0 ;
      end
%     WorkSpace: WK
%      wk = iniopt(wk,'a',     zeros(m1,n)) ;
      wk = iniopt(wk,'a',     0.0) ;
      if qrank1
          wk = iniopt(wk,'dxsave',zeros(n,nbroy)) ;
      else
          wk = iniopt(wk,'dxsave',0.0) ;
      end
      wk = iniopt(wk,'dx',    zeros(n,1)) ;
      wk = iniopt(wk,'dxq',   zeros(n,1)) ;
      wk = iniopt(wk,'xa',    zeros(n,1)) ;
      wk = iniopt(wk,'xwa',   zeros(n,1)) ;
      wk = iniopt(wk,'f',     zeros(n,1)) ;
      wk = iniopt(wk,'fa',    zeros(n,1)) ;
      wk = iniopt(wk,'eta',   zeros(n,1)) ;
      wk = iniopt(wk,'xw',    zeros(n,1)) ;
      wk = iniopt(wk,'fw',    zeros(n,1)) ;
      wk = iniopt(wk,'dxqa',  zeros(n,1)) ;
				
      wk = iniopt(wk,'sumxa0',0.0) ;
      wk = iniopt(wk,'sumxa1',0.0) ;
      wk = iniopt(wk,'fcmon',0.0) ;
      wk = iniopt(wk,'fcmin',0.0) ;
      wk = iniopt(wk,'sigma',0.0) ;
      wk = iniopt(wk,'sigma2',0.0) ;
      wk = iniopt(wk,'fca',0.0) ;
      wk = iniopt(wk,'fckeep',0.0) ;
      wk = iniopt(wk,'fcpri',0.0) ;
      wk = iniopt(wk,'dmycor',0.0) ;
      wk = iniopt(wk,'conv',0.0) ;
      wk = iniopt(wk,'sumx',0.0) ;
      wk = iniopt(wk,'sumxs',0.0) ;
      wk = iniopt(wk,'dlevf',0.0) ;
      wk = iniopt(wk,'niter',0) ;
      wk = iniopt(wk,'ncorr',0) ;
      wk = iniopt(wk,'nfcn',0) ;
      wk = iniopt(wk,'njac',0) ;
      wk = iniopt(wk,'nfcnj',0) ;
      wk = iniopt(wk,'nrejr1',0) ;
      wk = iniopt(wk,'new',0) ;
      wk = iniopt(wk,'iconv',0) ;
      
      iopt = iniopt(iopt,'norowscal',0) ;
	  if qinimo
		fprintf(lumon,'\n N = %4i\n\n Prescribed relative precision%10.2d\n',n,rtol) ;
		if jacgen == 1
		  jacg = 'a user function' ;
		elseif jacgen == 2
		  jacg = 'numerical differentiation (without feedback strategy)' ;
		elseif jacgen == 3
		  jacg = 'numerical differentiation (feedback strategy included)' ;
		end
		fprintf(lumon,'\n The Jacobian is supplied by %s\n',jacg) ;
		if mstor == 0
		  jacmode = 'full' ;
		elseif mstor == 1
		  jacmode = 'banded' ;
		end
		fprintf(lumon,' The Jacobian will be stored in %s mode\n',jacmode) ;
		if mstor == 1
		  fprintf(lumon,' Lower bandwidth : %3i   Upper bandwidth : %3i\n', ml,mu) ;
		end
		if iopt.norowscal == 1
		  rsmode = 'inhibited' ;
		else
		  rsmode = 'allowed' ;
		end
		fprintf(lumon,' Automatic row scaling of the jacobian is %s\n',rsmode) ;
	  end
	  [nonlin,iopt]=getopt(iopt,'nonlin',3) ;
	  iopt = iniopt(iopt,'boundeddamp',0) ;
	  if iopt.boundeddamp == 0
		qbdamp = nonlin == 4 ;
	  elseif iopt.boundeddamp == 1
		qbdamp = 1 ;
	  elseif iopt.boundeddamp == 2
		qbdamp = 0 ;
	  end
	  wk = iniopt(wk,'fcbnd',0.0) ;
	  if qbdamp
		if wk.fcbnd < 1.0
		   wk.fcbnd = 10.0 ;
		end
	  end
      if qinimo
        if qrank1
          rk1mode = 'allowed' ;
        else
          rk1mode = 'inhibited' ;
        end
        fprintf(lumon,'\n Rank-1 updates are %s\n',rk1mode) ;
        if nonlin == 1
          nonlinmode = 'linear' ;
        elseif nonlin == 2
          nonlinmode = 'mildly nonlinear' ;
        elseif nonlin == 3
          nonlinmode = 'highly nonlinear' ;
        elseif nonlin == 4
          nonlinmode = 'extremely nonlinear' ;
        end
        fprintf(lumon,' Problem is specified as being %s\n',nonlinmode) ;
        if qbdamp
          fprintf(lumon, ...
                  ' Bounded damping strategy is active\n bounding factor is %10.3e\n', wk.fcbnd) ;
        else
          fprintf(lumon, ' Bounded damping strategy is off\n') ;
        end
        if qordi 
          fprintf(lumon, ' Special mode: Ordinary Newton iteration will be done\n') ;
        end
        if qsimpl
           fprintf(lumon, ' Special mode: Simplified Newton iteration will be done\n') ;
        end
      end
%       Maximum permitted number of iteration steps
	  [nitmax,wk]=getopt(wk,'nitmax',50) ;
	  if nitmax <= 0
		nitmax=50 ;
	  end
	  wk.nitmax=nitmax ;
	  if qinimo
        fprintf(lumon,' Maximum permitted number of iteration steps : %6i\n',nitmax) ;
      end
%     Initial damping factor for highly nonlinear problems
	  wk = iniopt(wk,'fcstart',0.0) ;
	  qfcstr = wk.fcstart > 0.0 ;
	  if  ~ qfcstr
		wk.fcstart=1.0e-2 ;
		if nonlin == 4
		  wk.fcstart=1.0d-4 ;
		end
	  end
%     Minimal permitted damping factor
      wk = iniopt(wk,'fcmin',0.0) ;
	  if wk.fcmin <= 0.0
		wk.fcmin=1.0e-4 ;
		if nonlin == 4
		  wk.fcmin=1.0e-8 ;
		end
	  end
      fcmin=wk.fcmin ;
%     Rank1 decision parameter SIGMA
      wk = iniopt(wk,'sigma',0.0) ;
	  if wk.sigma < 1.0
	    wk.sigma=3.0 ;
	  end
	  if  ~ qrank1
	    wk.sigma=10.0/fcmin ;
	  end
%     decision parameter about increasing too small predictor
%     to greater corrector value
      wk = iniopt(wk,'sigma2',0.0) ;
	  if wk.sigma2 < 1.0
	    wk.sigma2=10.0/fcmin ;
	  end     
%     starting value of damping factor (fcmin <= fc <= 1.0)
	  if nonlin <= 2 &&  ~ qfcstr
%       for linear or mildly nonlinear problems
		fc = 1.0 ;
	  else
%       for highly or extremely nonlinear problems
		fc = wk.fcstart ;
	  end
%     simplied newton iteration implies ordinary newton it. mode
	  if iopt.qsimpl == 1
	    iopt.qordi = 1 ;
	  end
%     if ordinary newton iteration, damping factor is always one
	  if iopt.qordi == 1
	    fc = 1.0 ;
	  end
	  wk.fcstart = fc ;
	  if mprmon >= 2 &&  ~ qsucc
		fprintf(lumon,'\n\n%s\n\n%s%9.2e\n%s%9.2e\n%s%9.2e\n',' Internal parameters:', ...
		        ' Starting value for damping factor FCSTART = ', wk.fcstart, ...
		        ' Minimum allowed damping factor FCMIN = ', fcmin, ...
		        ' Rank-1 updates decision parameter SIGMA = ', wk.sigma ) ;
	  end
%
%       Initialize and start time measurements monitor
%
	  if  (~ qsucc)  &&  mprtim ~= 0 
	    dummy = 1 ;
		monini(' NLEQ1',lutim) ;
		mondef(0,'NLEQ1') ;
		mondef(1,'FCN     ') ;
		mondef(2,'Jacobi  ') ;
		mondef(3,'Lin-Fact') ;
		mondef(4,'Lin-Sol ') ;
		mondef(5,'Output  ') ;
		monsrt(dummy) ;
	  end
%
	  ierr=-1 ;
%       If IERR is unmodified on exit, successive steps are required
%       to complete the Newton iteration
	  if nbroy == 0
		 nbroy=1 ;
	  end
	  [x,xscal,rtol,iopt,ierr,xiter,sumxall,dlevfall,sumxqall,tolall,fcall,  ...
	            wk.a,wk.dxsave,wk.dx,wk.dxq,wk.xa,wk.xwa,wk.f,                                    ...
				wk.fa,wk.eta,wk.xw,wk.fw,wk.dxqa,wk.sumxa0,wk.sumxa1,wk.fcmon,wk.fcstart,         ...
				wk.fca,wk.fckeep,wk.fcpri,wk.dmycor,wk.conv,wk.sumx,wk.sumxs,wk.dlevf,            ...
				wk.niter,wk.ncorr,wk.nfcn,wk.njac,wk.nfcnj,wk.nrejr1,wk.new,wk.iconv,wk.qsucc] =  ...                            ...                    ...
		  n1int(n,fcn,par,x,xscal,rtol,nitmax,nonlin,iopt,ierr,wk,m1,m2,nbroy,                    ...
		        xiter,sumxall,dlevfall,sumxqall,tolall,fcall,   ...
				wk.a,wk.dxsave,wk.dx,wk.dxq,wk.xa,wk.xwa,wk.f,                                    ...
				wk.fa,wk.eta,wk.xw,wk.fw,wk.dxqa,wk.sumxa0,wk.sumxa1,wk.fcmon,wk.fcstart,         ...
				wk.fcmin,wk.sigma,wk.sigma2,wk.fca,wk.fckeep,                        ...
				wk.fcpri,wk.dmycor,wk.conv,wk.sumx,wk.sumxs,wk.dlevf,                ...
				mstor,mprerr,mprmon,mprsol,luerr,lumon,lusol,                        ...
				wk.niter,wk.ncorr,wk.nfcn,wk.njac,wk.nfcnj,wk.nrejr1,wk.new,wk.iconv,qbdamp) ;
%
	  if mprtim ~= 0 && ierr ~= -1 && ierr ~= 10
	    dummy = 1 ;
		monend(dummy) ;
	  end
%     Print statistics
      if mprmon >= 1 && ierr ~= -1 && ierr ~= 10
        fprintf(lumon,'\n%s%s%s\n%s%7i%s\n%s%7i%s\n%s%7i%s\n%s%7i%s\n%s%7i%s\n%s%7i%s\n%s\n\n', ...
                '   ******  Statistics * ',prodct, ' *******', ...
                '   ***  Newton iterations : ',wk.niter,  '  ***', ...
                '   ***  Corrector steps   : ',wk.ncorr,  '  ***', ...
                '   ***  Rejected rk-1 st. : ',wk.nrejr1, '  ***', ...
                '   ***  Jacobian eval.    : ',wk.njac,   '  ***',  ...
                '   ***  Function eval.    : ',wk.nfcn,   '  ***', ...
                '   ***  ...  for Jacobian : ',wk.nfcnj,  '  ***', ...
                '   *************************************') ;
      end
      info.xscal = xscal ;
      if ierr == -1
        info.rtol = tolall(wk.niter) ;
      else
        info.rtol = rtol ;
      end
      info.ierr = ierr ;
      info.xiter = xiter ;
      info.natlevel = sumxall ;
      info.simlevel = sumxqall ;
      info.stdlevel = dlevfall ;
      info.precision = tolall ;
      info.dampingfc = fcall ;
      info.niter = wk.niter ;
      info.ncorr = wk.ncorr ;
      info.nrejr1 = wk.nrejr1 ;
      info.njac = wk.njac ;
      info.nfcn = wk.nfcn ;
      info.nfcnj = wk.nfcnj ;
%     End of function NLEQ1

function [ierr,iopt,rtol,xscal] = n1pchk(n,x,xscal,rtol,iopt,wk)
%     ------------------------------------------------------------
%
%*    Summary :
%
%     N 1 P C H K : Checking of input parameters and options
%                   for NLEQ1.
%
%*    Parameters:
%     ===========
%
%     See parameter description in driver routine.
%
%*    functions called: D1MACH
%
%*    Machine dependent constants used:
%     =================================
%
%     EPMACH = relative machine precision
%     GREAT = squareroot of (maxreal divided by 10)
%     SMALL = squareroot of (smallest positive machine number
%             divided by relative machine precision)
%     ------------------------------------------------------------
      numopt=50 ;
%
%      ioptl = [0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, ...
%               0,0,0,0,0,0,0,0,0,0,-9999999,-9999999,-9999999,-9999999,-9999999      ] ;
%      ioptu = [1,1,3,1,0,9999999,9999999,0,1,0,3,99,6,99,3,99,0,0,1,     ...
%               99,0,0,0,0,0,0,0,0,0,0,4,1,1,1,1,0,0,2,3,0,0,0,0,0,0,     ...
%               9999999,9999999,9999999,9999999,9999999] ,
%
      epmach = d1mach(3) ;
      great  = sqrt(d1mach(2)/10.0) ;
      small  = d1mach(6) ;
      ierr = 0 ;
%        print error messages?
      mprerr = iopt.mprerr ;
      luerr  = iopt.luerr  ;
      if luerr  <=  0  ||  luerr  >  99
        luerr = 1 ;
        iopt.luerr=luerr ;
      end
%
%     checking dimensional parameter n
      if  n <= 0 
        if mprerr >= 1
          fprintf(luerr,'\n%s\n%s%5i\n',                                   ...
                  ' error: bad input to dimensional parameter n supplied', ...
                  '        choose n positive, your input is: n = ',n) ;
        end
        ierr = 20 ;
      end
%
%     problem type specification by user
      [nonlin,iopt] = getopt(iopt,'nonlin',0) ;
      if nonlin == 0
        nonlin = 3 ;
      end
      iopt.nonlin = nonlin ;
%
%     Checking and conditional adaption of the user-prescribed RTOL
      if rtol <= 0.0
        if mprerr >= 1
          fprintf(luerr,'\n%s\n',' Error: Nonpositive RTOL supplied') ;
        end
        ierr = 21 ;
      else
        tolmin = epmach*10.0*n ;
        if rtol < tolmin
          rtol = tolmin ;
          if mprerr >= 2
            fprintf(luerr, ...
            '\n warning: user prescribed rtol %s to reasonable %s value rtol = %11.2e\n', ...
            'increased','smallest',rtol ) ;
          end
        end
        tolmax = 1.0e-1 ;
        if rtol > tolmax
          rtol = tolmax ;
          if mprerr >= 2 
            fprintf(luerr, ...
            '\n warning: user prescribed rtol %s to reasonable %s value rtol = %11.2e\n', ...
            'decreased','largest',rtol) ;
          end
        end
      end
%     Test user prescribed accuracy and scaling on proper values
      if n <= 0
        return ;
      end
      if nonlin >= 3
        defscl = rtol ;
      else
        defscl = 1.0 ;
      end
      for i=1:n
        if xscal(i) < 0.0
          if mprerr >= 1
            fprintf(luerr,' error: negative value in xscal(%5i) supplied\n', i) ;
          end
          ierr = 22 ;
        end
        if xscal(i) == 0.0
          xscal(i) = defscl ;
        end
        if  xscal(i) > 0.0  &&  xscal(i) < small 
          if mprerr >= 2
            fprintf(luerr,' warning: xscal(%5i) = %9.2e too small, increased to %9.2e\n', ...
                    i,xscal(i),small ) ;
          end
          xscal(i) = small ;
        end
        if xscal(i) > great
          if mprerr >= 2
            fprintf(luerr,' warning: xscal(%5i) = %9.2e too big, decreased to %9.2e\n', ...
                    i,xscal(i),great ) ;
          end
          xscal(i) = great ;
        end
      end
%     Special dependence on Jacobian's storage mode for ML and MU
%      [mstor,iopt] = getopt(iopt,'mstor',0) ;
%      if mstor == 0
%        ioptu(6)=0 ;
%        ioptu(7)=0 ;
%      elseif mstor == 1
%        ioptu(6)=n-1 ;
%        ioptu(7)=n-1 ;
%      end
%     Checks options
%      DO 20 I=1,30
%        if IOPT(I) < IOPTL(I)  ||  IOPT(I) > IOPTU(I)
%          IERR=30
%          if MPRERR >= 1
%            WRITE(LUERR,20001) I,IOPT(I),IOPTL(I),IOPTU(I)
%20001       FORMAT(' Invalid option specified: IOPT(',I2,')=',I12,';',
%     $             /,3X,'range of permitted values is ',I8,' to ',I8)
%          end
%        end
%20    CONTINUE
%     End of function N1PCHK

function [x,xscal,rtol,iopt,ierr,xiter,sumxall,dlevfall,sumxqall,tolall,fcall,         ...
          a,dxsave,dx,dxq,xa,xwa,f,fa,eta,xw,fw,dxqa,                                  ...
          sumxa0,sumxa1,fcmon,fc,fca,fckeep,fcpri,dmycor,conv,sumx,sumxs,dlevf,        ...
          niter,ncorr,nfcn,njac,nfcnj,nrejr1,new,iconv,qsucc] =                        ...
     n1int(n,fcn,par,x,xscal,rtol,nitmax,nonlin,iopt,ierr,wk,m1,m2,nbroy,         ...
           xiter,sumxall,dlevfall,sumxqall,tolall,fcall,        ...
           a,dxsave,dx,dxq,xa,xwa,f,fa,eta,xw,fw,dxqa,sumxa0,sumxa1,fcmon,        ...
           fc,fcmin,sigma,sigma2,fca,fckeep,fcpri,dmycor,conv,sumx,sumxs,dlevf,   ...
           mstor,mprerr,mprmon,mprsol,luerr,lumon,lusol,                          ...
           niter,ncorr,nfcn,njac,nfcnj,nrejr1,new,iconv,qbdamp)
%     ------------------------------------------------------------
%
%*    Summary :
%
%     N 1 I N T : Core function for NLEQ1 .
%     Damped Newton-algorithm for systems of highly nonlinear
%     equations especially designed for numerically sensitive
%     problems.
%
%*    Parameters:
%     ===========
%
%       N,FCN,X,XSCAL,RTOL   
%                         See parameter description in driver routine
%
%       NITMAX      Int    Maximum number of allowed iterations
%       NONLIN      Int    Problem type specification
%                          (see IOPT-field NONLIN)
%       IOPT        Int    See parameter description in driver routine
%       IERR        Int    See parameter description in driver routine
%       WK          Struct Workspace structure array
%       M1          Int    Leading dimension of Jacobian array A
%                          for full case Jacobian: N
%                          for banded Jacobian: 2*ML+MU+1;
%                          ML, MU see IOPT-description in driver routine
%       M2          Int    for full case Jacobian: N
%                          for banded Jacobian: ML+MU+1
%       NBROY       Int    Maximum number of possible consecutive
%                          iterative Broyden steps. (See wk.nbroy)
%       A(M1,N)     Dble   Holds the Jacobian matrix.
%       DXSAVE(X,NBROY)
%                   Dble   Used to save the quasi Newton corrections of
%                          all previously done consecutive Broyden
%                          steps.
%       DX(N)       Dble   Current Newton correction
%       DXQ(N)      Dble   Simplified Newton correction J(k-1)*X(k)
%       XA(N)       Dble   Previous Newton iterate
%       XWA(N)      Dble   Scaling factors used for latest decomposed
%                          Jacobian for column scaling - may differ
%                          from XW, if Broyden updates are performed
%       F(N)        Dble   Function (FCN) value of current iterate
%       FA(N)       Dble   Function (FCN) value of previous iterate
%       ETA(N)      Dble   Jacobian approximation: updated scaled
%                          denominators
%       XW(N)       Dble   Scaling factors for iteration vector
%       FW(N)       Dble   Scaling factors for rows of the system
%       DXQA(N)     Dble   Previous Newton correction
%       FC          Dble   Current Newton iteration damping factor.
%       FCMIN       Dble   Minimum permitted damping factor. If
%                          FC becomes smaller than this value, one
%                          of the following may occur:
%                          a.    Recomputation of the Jacobian
%                                matrix by means of difference
%                                approximation (instead of Rank1
%                                update), if Rank1 - update
%                                previously was used
%                          b.    Fail exit otherwise
%       SIGMA       Dble   Decision parameter for rank1-updates
%       SIGMA2      Dble   Decision parameter for damping factor
%                          increasing to corrector value
%       FCA         Dble   Previous Newton iteration damping factor.
%       FCKEEP      Dble   Keeps the damping factor as it is at start
%                          of iteration step.
%       CONV        Dble   Scaled maximum norm of the Newton-
%                          correction. Passed to RWK-field on output.
%       SUMX        Dble   Square of the natural level (see equal-
%                          named RWK-output field)
%       SUMXS       Dble   Square of the "simplified" natural level
%                          (see equal-named RWK-internal field)
%       DLEVF       Dble   The standard level (see equal-
%                          named RWK-output field)
%       MSTOR       Int    see description of IOPT-field MSTOR
%       MPRERR,MPRMON,MPRSOL,LUERR,LUMON,LUSOL :
%                          See description of equal named IOPT-fields
%       NITER,NCORR,NFCN,NJAC,NFCNJ,NREJR1,NEW :
%                          See description of equal named WK-fields
%                          in the driver function
%       QBDAMP      Logic  Flag, that indicates, whether bounded damping
%                          strategy is active:
%                          1 = bounded damping strategy is active
%                          0 = normal damping strategy is active
%
%*    Internal double variables
%     =========================
%
%       AJDEL    See wk.ajdel (num. diff. without feedback)
%       AJMIN    See wk.ajmin (num. diff. without feedback)
%       CONVA    Holds the previous value of CONV .
%       DMUE     Temporary value used during computation of damping 
%                factors predictor.
%       EPDIFF   sqrt(10*epmach) (num. diff. with feedback)
%       ETADIF   See description of wk.etadif (num. diff. with feedback)
%       ETAINI   Initial value for all ETA-components (num. diff. fb.)
%       ETAMAX   Maximum allowed pertubation (num. diff. with feedback)
%       ETAMIN   Minimum allowed pertubation (num. diff. with feedback)
%       FCDNM    Used to compute the denominator of the damping 
%                factor FC during computation of it's predictor,
%                corrector and aposteriori estimate (in the case of
%                performing a Rank1 update) .
%       FCK2     Aposteriori estimate of FC.
%       FCH      Temporary storage of new corrector FC.
%       FCMIN2   FCMIN^2 . Used for FC-predictor computation.
%       FCNUMP   Gets the numerator of the predictor formula for FC.
%       FCNMP2   Temporary used for predictor numerator computation.
%       FCNUMK   Gets the numerator of the corrector computation 
%                of FC .
%       SUMXA    Natural level of the previous iterate.
%       TH       Temporary variable used during corrector- and 
%                aposteriori computations of FC.
%
%*    Internal integer variables
%     ==========================
%
%     IFAIL      Gets the return value from functions called from
%                N1INT (N1FACT, N1SOLV, FCN) 
%     ISCAL      Holds the scaling option from iopt.iscal      
%     MODE       Matrix storage mode (see iopt.mode) 
%     NRED       Count of successive corrector steps
%
%
%*    Internal logical variables
%     ==========================
%
%     QGENJ      Jacobian updating technique flag:
%                =1 : Call of analytical function FCN with Jacobian
%                     flag or numerical differentiation
%                =0 : rank1- (Broyden-) update
%     QINISC     Iterate initial-scaling flag:
%                =1 : at first call of N1SCAL
%                =0 : at successive calls of N1SCAL
%     QSUCC      See description of IOPT-field QSUCC.
%     QJCRFR     Jacobian refresh flag:
%                set to 1 if damping factor gets too small
%                and Jacobian was computed by rank1-update. 
%                Indicates, that the Jacobian needs to be recomputed
%                by user function FCN or numerical differentation.
%     QLINIT     Initialization state of linear solver workspace:
%                =0 : Not yet initialized
%                =1 : Initialized - N1FACT has been called at
%                     least one time.
%     QSCALE     Holds the value of  ~ norowscal. See description
%                of IOPT-field norowscal.
%
%*    Functions called:
%     =================
%
%       N1FACT, N1SOLV, N1JAC,  N1JACB, N1JCF,  N1JCFB, N1LVLS, 
%       N1SCRF, N1SCRB, N1SOUT, N1PRV1, N1PRV2, N1SCAL,
%       MONON,  MONOFF, D1MACH, WNORM
%
%     ------------------------------------------------------------
%
      persistent clin0 clin1 calpha alphae alphak alphaa sumxa2 qmstop l u p;
      one=1.0 ; zero=0.0 ; half=0.5 ; ten=10.0 ;
      epmach = d1mach(3) ;
      small  = d1mach(6) ;
      n = length(x) ;
%*    Begin
%       ----------------------------------------------------------
%       1 Initialization
%       ----------------------------------------------------------
%       1.1 Control-flags and -integers
        qsucc = wk.qsucc ;
        qscale = iopt.norowscal ~= 1 ;
        qordi  = iopt.qordi ;
        qsimpl = iopt.qsimpl ;
        qrank1 = iopt.qrank1 ;
        [iormon,iopt] = getopt(iopt,'iormon',0) ;
        if iormon == 0 
          iormon=2 ;
        end
        [iscal,iopt] = getopt(iopt,'iscal',0) ;
        [mode,iopt]  = getopt(iopt,'mode' ,0) ;
        jacgen = iopt.jacgen ;
        qmixio = lumon == lusol && mprmon ~= 0 && mprsol ~= 0 ;
        qlu    = ~ qsimpl ;
        mprtim = iopt.mprtim ;
%       ----------------------------------------------------------
%       1.2 Derivated dimensional parameters
        if mstor == 0
          ml=0 ;
          mu=0 ;
        elseif mstor == 1
          ml=m1-m2 ;
          mu=m2-1-ml ;
        end
%       ----------------------------------------------------------
%       1.3 Derivated internal parameters
        fcmin2 = fcmin*fcmin ;
        rsmall = sqrt(ten*rtol) ;
%       ----------------------------------------------------------
%       1.4 Adaption of input parameters, if necessary
        if fc < fcmin
          fc = fcmin ;
        end
        if fc > one
          fc = one ;
        end
%       ----------------------------------------------------------
%       1.5 Initial preparations
        qjcrfr = 0 ;
        qlinit = 0 ;
        qiter = 1 ;
        ifail = 0 ;
        fcbnd = zero ;
        if qbdamp
          fcbnd = wk.fcbnd ;
        end
%       ----------------------------------------------------------
%       1.5.1 Numerical differentiation related initializations
        if jacgen == 2
          [ajdel,wk] = getopt(wk,'ajdel',0.0) ;
          if ajdel <= small
            ajdel = sqrt(epmach*ten) ;
          end
          [ajmin,wk] = getopt(wk,'ajmin',0.0) ;
        elseif jacgen == 3
          [etadif,wk] =  getopt(wk,'etadif',0.0) ;
          if etadif  <=  small
            etadif = 1.0e-6 ;
          end
          [etaini,wk] =  getopt(wk,'etaini',0.0) ;
          if etaini  <=  small
            etaini = 1.0e-6 ;
          end
          epdiff = sqrt(epmach*ten) ;
          etamax = sqrt(epdiff) ;
          etamin = epdiff*etamax ;
        end
%       ----------------------------------------------------------
%       1.5.2 Miscellaneous preparations of first iteration step
        if  ~ qsucc
          niter = 0 ;
          wk.niter = niter ;
          xiter(1:n,1) = x ;
          ncorr = 0 ;
          nrejr1 = 0 ;
          nfcn = 0 ;
          njac = 0 ;
          nfcnj = 0 ;
          qgenj = 1 ;
          qinisc = 1 ;
          fckeep = fc ;
          fca = fc ;
          fcpri = fc ;
          fck2 = fc ;
          fcmon = fc ;
          conv = zero ;
          if jacgen == 3
            eta = etaini*ones(n,1) ;
          end
          xa = x ;
          iconv = 0 ;
          alphae = zero ;
          calpha = zero ;
          sumxa1 = zero ;
          sumxa0 = zero ;
          clin0  = zero ;
          qmstop = 0 ;
%         ------------------------------------------------------
%         1.6 Print monitor header
          if mprmon >= 2  && ~ qmixio
            fprintf(lumon,'\n\n\n%1s',' ') ;
            fprintf(lumon,'%s\n\n%s\n', ...
                  '  ******************************************************************' , ...
                  '        It       Normf           Normx         Damp.Fct.   New'        ) ;
          end
%         --------------------------------------------------------
%         1.7 Startup step
%         --------------------------------------------------------
%         1.7.1 Computation of the residual vector
          if mprtim ~= 0
             monon(1) ;
          end
          [f,ifail] = feval(fcn,x,'',par) ;
          if mprtim ~= 0
            monoff(1) ;
          end
          nfcn = nfcn+1 ;
%     Exit, if ...
          if ifail ~= 0
            ierr = 82 ;
            qiter = 0 ;
          end
          if length(f) ~= n
            ierr = 22 ;
            qiter = 0 ;
          end
        else
          qinisc = 0 ;
        end
%
%       Main iteration loop
%       ===================
%
%       Repeat
        while qiter
%         --------------------------------------------------------
%         2 Startup of iteration step
          if  ~ qjcrfr
%           ------------------------------------------------------
%           2.1 Scaling of variables X(N)
            xw = n1scal(n,x,xa,xscal,iscal,qinisc,iopt) ;
            qinisc = 0 ;
            if niter ~= 0
%             ----------------------------------------------------
%             2.2 Aposteriori estimate of damping factor
              dxqa = dxq ;
              if  ~ qordi
                fcnump = sum((dx./xw).^2) ;
                th = fc-one ;
                fcdnm = sum(((dxqa+th*dx)./xw).^2) ;
%               --------------------------------------------------
%               2.2.2 Decision criterion for Jacobian updating
%                     technique:
%                     QGENJ == 1 numerical differentation,
%                     QGENJ == 0 rank1 updating
                qgenj = 1 ;
                if fc == fcpri
                  qgenj = fc < one || fca < one || dmycor <= fc*sigma ...
                          || ~ qrank1 || new+2 > nbroy ;
                  fca = fc ;
                else
                  dmycor = fca*fca*half*sqrt(fcnump/fcdnm) ;
                  if nonlin <= 3
                    fccor = min(one,dmycor) ;
                  else
                    fccor = min(one,half*dmycor) ;
                  end
                  fca = max(min(fc,fccor),fcmin) ;
%$test-begin
                  if mprmon >= 5
                    fprintf(lumon,'%s\n%s%18.10e  %s%18.10e\n%s%18.10e  %s%18.10e\n%s%18.10e\n%s\n', ...
                      ' +++ aposteriori estimate +++', ...
                      '  fccor  = ', fccor,  ...
                      '  fc     = ', fc,     ...
                      '  dmycor = ', dmycor, ... 
                      '  fcnump = ', fcnump, ...
                      '  fcdnm  = ', fcdnm,  ...
                      ' ++++++++++++++++++++++++++++' ) ;
                  end
%$test-end 
                end
                fck2 = fca ;
%               ------------------------------------------------------
%               2.2.1 Computation of the numerator of damping
%                     factor predictor
                fcnmp2 = sum((dxqa./xw).^2) ;
                fcnump = fcnump*fcnmp2 ;
              end
            end
          end
          qjcrfr = 0 ;
%         --------------------------------------------------------
%         2.3 Jacobian matrix (stored to array A(M1,N))
%         --------------------------------------------------------
%         2.3.1 Jacobian generation by routine JAC or
%               difference approximation (If QGENJ == 1)
%               - or -
%               Rank-1 update of Jacobian (If QGENJ == 0)
          if qgenj  &&  ( ~ qsimpl  ||  niter == 0 )
            new = 0 ;
            if jacgen == 1
               if mprtim ~= 0
                 monon(2) ;
               end
               [a,ifail] = feval(fcn,x,'jacobian',par) ;
               if mprtim ~= 0
                 monoff(2) ;
               end
            else
              if mstor == 0
                if mprtim ~= 0
                  monon(2) ;
                end
                if jacgen == 3
                  [a,eta,nfcnj,ifail] = n1jcf(fcn,par,n,n,x,f,xw,eta,etamin,etamax,etadif,conv,nfcnj) ;
                end
                if jacgen == 2 
                  [a,nfcnj,ifail] = n1jac(fcn, par, n, n, x, f, xw, ajdel, ajmin, nfcnj) ;
                end
                if mprtim ~= 0
                  monoff(2) ;
                end
              elseif mstor == 1
                if mprtim ~= 0
                  monon(2) ;
                end
                if jacgen == 3 
                  [a,eta,nfcnj,ifail] = n1jcfb(fcn,par,n,m1,ml,x,f,xw,eta,etamin,etamax, ...
                                               etadif,conv,nfcnj) ;
                end
                if jacgen == 2 
                  [a,nfcnj,ifail] = n1jacb (fcn, par, n, m1, ml, x, f, xw, ajdel, ajmin, nfcnj) ;
                end
                if mprtim ~= 0
                  monoff(2) ;
                end
              end
            end
            njac = njac + 1 ;
%     Exit, If ...
            if jacgen == 1  &&  ifail < 0
              ierr = 83;
              break ;
            end
            if jacgen ~= 1  &&  ifail ~= 0
              ierr = 82 ;
              break ;
            end
          elseif  ~ qsimpl
            new = new+1 ;
          end
          if  new == 0  &&  (qlu || niter == 0) 
%             ------------------------------------------------------
%             2.3.2.1 Save scaling values
              xwa = xw(1:n) ;
%             ------------------------------------------------------
              if issparse(a)
                  nza = nnz(a) ;
                  [row,col] = find(a ~= 0) ;
                  for k=1:nza
                      a(row(k),col(k))=-a(row(k),col(k))*xw(col(k)) ;
                  end
              else
%             2.3.2.2 Prepare Jac. for use by band-solver DGBFA/DGBSL
                  if mstor == 1
                    for l1=1:n
                      for l2=m2:-1:1
                        a(l2+ml,l1)=a(l2,l1) ;
                      end
                    end
                  end
%             ------------------------------------------------------
%             2.4 Prepare solution of the linear system
%             ------------------------------------------------------
%             2.4.1 internal column scaling of matrix A
                  if mstor == 0
                    for k=1:n
                      a(1:n,k)=-a(1:n,k)*xw(k) ;
                    end
                  elseif mstor == 1
                    for k=1:n
                      l2=max(1+m2-k,ml+1) ;
                      l3=min(n+m2-k,m1)   ;
                      a(l2:l3,k) = -a(l2:l3,k)*xw(k) ;
                    end
                  end
              end
%             ----------------------------------------------------
%             2.4.2 Row scaling of matrix A
              if qscale
                if mstor == 0
                  [a,fw] = n1scrf(n,n,a) ;
                elseif mstor == 1
                  [a,fw] = n1scrb(n,m1,ml,mu,a) ;
                end
              else
                fw = ones(n,1) ;
              end
          end
%         --------------------------------------------------------
%         2.4.3 Save and scale values of F(N)
          fa = f ;
          t1 = f .* fw ;
%         --------------------------------------------------------
%         3 Central part of iteration step
%         --------------------------------------------------------
%         3.1 Solution of the linear system
%         --------------------------------------------------------
%         3.1.1 Decomposition of (N,N)-matrix A
          if  new == 0  &&  (qlu || niter == 0)
            if mprtim ~= 0
              monon(3) ;
            end
            [l,u,p,ifail] = n1fact(n,m1,ml,mu,a,iopt) ;
            if mprtim ~= 0
               monoff(3) ;
            end
%       exit repeat if ...
            if ifail ~= 0
              if ifail == 1
                ierr = 1 ;
              else
                ierr = 80 ;
              end
              break ;
            end
          end
          qlinit = 1 ;
%         --------------------------------------------------------
%         3.1.2 Solution of linear (N,N)-system
          if new == 0
            if mprtim ~= 0
              monon(4) ;
            end
            [t1,ifail] = n1solv(n,m1,ml,mu,l,u,p,t1,iopt) ;
            if mprtim ~= 0
              monoff(4) ;
            end
%     Exit Repeat If ...
            if ifail ~= 0
              ierr = 81 ;
              break ;
            end
          else  
            alfa1 = sum(dx.*dxq./xw.^2) ;
            alfa2 = sum(dx.^2./xw.^2) ;
            alfa=alfa1/alfa2 ;
            beta=one-alfa ;
            t1 = (dxq+(fca-one)*alfa*dx)/beta ;
            if new == 1
              dxsave(1:n,1) = dx ;
            end
            dxsave(1:n,new+1)=t1 ;
            dx=t1 ;
            t1=t1./xw ;
          end
%         --------------------------------------------------------
%         3.2 Evaluation of scaled natural level function SUMX
%             scaled maximum error norm CONV
%             evaluation of (scaled) standard level function
%             DLEVF ( DLEVF only, if MPRMON >= 2 )
%             and computation of ordinary Newton corrections 
%             DX(N)
          if  ~  qsimpl
            [dx,conv,sumx,dlevf] = n1lvls(n,dx,t1,xw,f,mprmon,new == 0) ;
          else
            [dx,conv,sumx,dlevf] = n1lvls(n,dx,t1,xwa,f,mprmon,new == 0) ;
          end
          wk.sumx = sumx ;  wk.dlevf = dlevf ;
          xa=x ;
          sumxa = sumx ;
          dlevxa = sqrt(sumxa/n) ;
          conva = conv ;
          dxanrm = wnorm(n,dx,xw) ;
          sumxall(niter+1) = dlevxa ;  dlevfall(niter+1) = dlevf ;
%         --------------------------------------------------------
%         3.3 A - priori estimate of damping factor FC
          if niter ~= 0 && nonlin ~= 1 && new == 0 && ~ qordi
%           ------------------------------------------------------
%           3.3.1 Computation of the denominator of a-priori
%                 estimate
            fcdnm = sum(((dx-dxqa)./xw).^2)*sumx ;
%           ------------------------------------------------------
%           3.3.2 New damping factor
            if fcdnm > fcnump*fcmin2  || (nonlin == 4  &&  fca^2*fcnump  <  4.0*fcdnm)
              dmypri = fca*sqrt(fcnump/fcdnm) ;
              fcpri = min(dmypri,one) ;
              if nonlin == 4
                fcpri = min(half*dmypri,one) ;
              end
            else
              fcpri = one ;
%$test-begin
              dmypri = -1.0 ;
%$Test-end
            end
%$Test-begin
            if mprmon >= 5
              fprintf(lumon,'%s\n%s%18.10e  %s%18.10e\n%s%18.10e  %s%18.10e\n%s%18.10e  %s%18.10e\n%s\n', ...
                      ' +++ apriori estimate +++', ...
                      '  fcpri  = ', fcpri,  ...
                      '  fc     = ', fc,     ...
                      '  fca    = ', fca,    ... 
                      '  dmypri = ', dmypri, ...
                      '  fcnump = ', fcnump, ...
                      '  fcdnm  = ', fcdnm,  ...
                      ' ++++++++++++++++++++++++++++' ) ;
            end
%$Test-end 
            fc = max(fcpri,fcmin) ;
            if qbdamp
              fcbh = fca*fcbnd ;
              if fc > fcbh
                fc = fcbh ;
                if mprmon >= 4
                  fprintf(lumon,'%s\n',' *** incr. rest. act. (a prio) ***') ;
                end
              end
              fcbh = fca/fcbnd ;
              if fc < fcbh
                fc = fcbh ;
                if mprmon >= 4
                  fprintf(lumon,'%s\n',' *** decr. rest. act. (a prio) ***'); 
                end
              end
            end
          end
%
          if iormon >= 2
            sumxa2 = sumxa1 ;
            sumxa1 = sumxa0 ;
            sumxa0 = dlevxa ;
            if sumxa0 == zero
              sumxa0 = small ;
            end ;
%           Check convergence rates (linear and superlinear)
%           ICONV : Convergence indicator
%                   =0: No convergence indicated yet
%                   =1: Damping factor is 1.0d0
%                   =2: Superlinear convergence detected (alpha >=1.2)
%                   =3: Quadratic convergence detected (alpha > 1.8)
            fcmon = min(fc,fcmon) ;
            if fcmon < one
              iconv = 0 ;
              alphae = zero ;
            end
            if fcmon == one  &&  iconv == 0
              iconv=1 ;
            end
            if niter >= 1
              clin1 = clin0 ;
              clin0 = sumxa0/sumxa1 ;
            end
            if iconv >= 1 && niter >= 2
              alphak = alphae ;
              alphae = zero ;
              if clin1 <= 0.95
                alphae = log(clin0)/log(clin1) ;
              end
              if alphak ~= zero
                alphak =0.5*(alphae+alphak) ;
              end
              alphaa = min(alphak,alphae) ;
              calphk = calpha ;
              calpha = zero ;
              if alphae ~= zero
                calpha = sumxa1/sumxa2^alphae ;
              end
              sumxte = sqrt(calpha*calphk)*sumxa1^alphak-sumxa0 ;
              if alphaa >= 1.2  &&  iconv == 1
                iconv = 2 ;
              end
              if alphaa > 1.8
                iconv = 3 ;
              end
              if mprmon >= 4
                fprintf(lumon,'%s%1i  %s%9.2e  %s%9.2e  %s%9.2e%s\n%s%9.2e  %s%9.2e%s\n', ...
                         ' ** iconv: ',                iconv, ...
                         '  alpha: ',                  alphae, ...
                         '  const-alpha: ',            calpha, ...
                         '  const-lin: ',              clin0,    ' **', ...
                         ' **           alpha-post: ', alphak, ...
                         ' check: ',                   sumxte, ...
                         '                         **')
              end
              if  iconv >= 2  &&  alphaa < 0.9
                 if iormon == 3
                   ierr = 4 ;
                   break ;
                 else
                   qmstop = 1 ;
                 end 
              end
            end
          end
          fcmon = fc ;
%
%         --------------------------------------------------------
%         3.4 Save natural level for later computations of
%             corrector and print iterate
          fcnumk = sumx ;
          if mprmon >= 2
            if mprtim ~= 0
              monon(5) ;
            end
            n1prv1(dlevf,dlevxa,fckeep,niter,new,mprmon,lumon,qmixio) ;
            if mprtim ~= 0
              monoff(5) ;
            end
          end
          nred = 0 ;
          qnext = 0 ;
          qrep  = 0 ;
          qred  = 1 ;
          icnv = 0 ;
%
%         Damping-factor reduction loop
%         ================================
          while qred
%           ------------------------------------------------------
%           3.5 Preliminary new iterate
            x=xa+dx*fc ;
            fcall(niter+1) = fc ;
%           -----------------------------------------------------
%           3.5.2 Exit, if problem is specified as being linear
%     Exit Repeat If ...
            if  nonlin == 1 
              ierr = 0 ;
              break ;
            end
%           ------------------------------------------------------
%           3.6.1 Computation of the residual vector
            if mprtim ~= 0
              monon(1) ;
            end
            [f,ifail] = feval(fcn,x,'',par) ;
            if mprtim ~= 0
              monoff(1) ;
            end
            nfcn = nfcn+1 ;
%     exit, if ...
            if ifail < 0
              ierr = 82 ;
              break ; 
            end
            if ifail == 1  ||  ifail == 2
              if qordi
                ierr = 82 ;
                break ;
              end
              if ifail == 1
                fcredu = half ;
              else
                fcredu = f(1) ;
%     exit, if ...
                if fcredu <= 0  ||  fcredu >= 1
                  ierr = 82 ;
                  break ;
                end
              end
              if mprmon >= 2
                fprint(lumon,'%8s%2i%37s%7.5f%4s%2i\n',               ...
                             '        ',niter,                        ...
                             ' fcn could not be evaluated          ', ...
                             fc,'    ',new ) ;
              end
              fch = fc ;
              fc = fcredu*fc ;
              if fch > fcmin
                fc = max(fc,fcmin) ;
              end
              if qbdamp
                fcbh = fch/fcbnd ;
                if fc < fcbh
                  fc = fcbh ;
                  if mprmon >= 4
                    fprint(lumon,'%s\n',' *** decr. rest. act. (fcn redu.) ***') ;
                  end
                end
              end
              if fc < fcmin
                ierr = 3 ;
                break ;
              end  
            else
              if qordi
%             -------------------------------------------------------
%             3.6.2 Convergence test for ordinary Newton iteration
%     Exit Repeat If ...
                tolall(niter+1) = dxanrm ;
                if  dxanrm <= rtol 
                  ierr = 0 ;
                  icnv = 1 ;
                  break ;
                end
              else
                t1 = f .* fw ;
%               --------------------------------------------------
%               3.6.3 Solution of linear (N,N)-system
                if mprtim ~= 0
                  monon(4) ;
                end
                [t1,ifail] = n1solv(n,m1,ml,mu,l,u,p,t1,iopt) ;
                if mprtim ~= 0
                  monoff(4) ;
                end
%     Exit Repeat If ...
                if ifail ~= 0
                  ierr = 81 ;
                  break ;
                end
                if new > 0
                  dxq = t1.*xwa ;
                  for iloop=1:new 
                    sum1 = sum((dxq.*dxsave(1:n,iloop))./ xw.^2) ;
                    sum2 = sum((dxsave(1:n,iloop)./xw).^2) ;
                    beta = sum1/sum2 ;
                    dxq = dxq+beta*dxsave(1:n,iloop+1) ;
                    t1 = dxq./xw ;
                  end
                end
%               ----------------------------------------------------
%               3.6.4 Evaluation of scaled natural level function
%                     SUMX
%                     scaled maximum error norm CONV and evaluation
%                     of (scaled) standard level function DLEVFN
                if  ~  qsimpl
                  [dxq,conv,sumx,dlevfn] = n1lvls(n,dxq,t1,xw,f,mprmon,new == 0) ;
                else
                  [dxq,conv,sumx,dlevfn] = n1lvls(n,dxq,t1,xwa,f,mprmon,new == 0) ;
                end 
                sumxqall(niter+1) = sqrt(sumx/n) ; 
                dxnrm = wnorm(n,dxq,xw) ;
%               -----------------------------------------------------
%               3.6.5 Convergence test
%     Exit Repeat If ...
                tolall(niter+1) = dxnrm ;
                if  dxnrm <= rtol  &&  dxanrm <= rsmall  &&  fc == one 
                  ierr = 0 ;
                  icnv = 1 ;
                  break ;
                end
%           
                fca = fc ;
%               ----------------------------------------------------
%               3.6.6 Evaluation of reduced damping factor
                th = fca-one ;
                fcdnm=sum(((dxq+th*dx)./xw).^2) ;
                if fcdnm ~= zero
                  dmycor = fca*fca*half*sqrt(fcnumk/fcdnm) ;
                else
                  dmycor = 1.0e+35 ;
                end
                if nonlin <= 3
                  fccor = min(one,dmycor) ;
                else
                  fccor = min(one,half*dmycor) ;
                end
%$Test-begin
                if mprmon >= 5
                  fprintf(lumon, ...
                      '%s\n%s%18.10e  %s%18.10e\n%s%18.10e  %s%18.10e\n%s%18.10e  %s%18.10e\n%s\n', ...
                      ' +++ corrector computation +++', ...
                      '  fccor  = ', fccor,  ...
                      '  fc     = ', fc,     ...
                      '  dmycor = ', dmycor, ...
                      '  fcnumk = ', fcnumk, ...
                      '  fcdnm  = ', fcdnm,  ...
                      '  fca    = ', fca,    ... 
                      ' ++++++++++++++++++++++++++++' ) ;
                end
%$Test-end 
              end
%             ------------------------------------------------------
%             3.7 Natural monotonicity test
              if sumx > sumxa  && ~ qordi
%               ----------------------------------------------------
%               3.8 Output of iterate
                if mprmon >= 3
                  if mprtim ~= 0
                    monon(5) ;
                  end
                  n1prv2(dlevfn,sqrt(sumx/n),fc,niter,mprmon,lumon,qmixio,'*') ;
                  if mprtim ~= 0
                    monoff(5) ;
                  end
                end
                if qmstop
                  ierr = 4 ;
                  break ;
                end
                fch = min(fccor,half*fc) ;
                if fc > fcmin
                  fc=max(fch,fcmin) ;
                else
                  fc=fch ;
                end
                if qbdamp
                  fcbh = fca/fcbnd ;
                  if fc < fcbh
                    fc = fcbh ;
                    if mprmon >= 4
                      fprintf(lumon,'%s\n',' *** decr. rest. act. (a post) ***') ;
                    end
                  end
                end
                fcmon = fc ;
%  
%$Test-begin
                  if mprmon >= 5
                    fprintf(lumon,'%s\n%s%18.10e\n%s\n',    ...
                      ' +++ corrector setting 1 +++',   ...
                      ' fc     = ', fc,                 ...
                      ' +++++++++++++++++++++++++++') ;
                  end
%$Test-end 
                qrep = 1 ;
                ncorr = ncorr+1 ;
                nred = nred+1 ;
%               ----------------------------------------------------
%               3.10 If damping factor is too small:
%                    Refresh Jacobian,if current Jacobian was computed
%                    by a Rank1-update, else fail exit
                qjcrfr  = fc < fcmin || (new > 0 && nred > 1) ;
%     Exit Repeat If ...
                if qjcrfr && new == 0
                  ierr = 3 ;
                  break ;
                end
              else
                if  ~ qordi && ~ qrep  &&  fccor > sigma2*fc
                  if mprmon >= 3
                    if mprtim ~= 0
                      monon(5) ;
                    end
                    n1prv2(dlevfn,sqrt(sumx/dble(float(n))),fc, ...
                            niter,mprmon,lumon,qmixio,'+')
                    if mprtim ~= 0
                      monoff(5) ;
                    end
                  end
                  fc = fccor ;
%$test-begin
                  if mprmon >= 5
                    fprintf(lumon,'%s\n%s%18.10e\n%s\n',    ...
                      ' +++ corrector setting 2 +++',   ...
                      ' fc     = ', fc,                 ...
                      ' +++++++++++++++++++++++++++') ;
                  end
%$test-end 
                  qrep = 1 ;
                else
                  qnext = 1 ;
                end
              end
%3109     continue
            end
            qred =  ~ (qnext || qjcrfr) ;
          end
%         UNTIL ( expression - negated above)
%         End of damping-factor reduction loop
          if nonlin == 1 || icnv == 1 || (ierr ~= 0 && ierr ~= -1)
            break ;
          end
%         =======================================
          if qjcrfr
%           ------------------------------------------------------
%           3.11 Restore former values for repeting iteration
%                step
            nrejr1 = nrejr1+1 ;
            x = xa ;
            f = fa ;
            if mprmon >= 2
              fprint(lumon,'%8s%2i%37s%7.5f%4s%2i\n',                 ...
                           '        ',niter,                        ...
                           ' not accepted damping factor         ', ...
                           fc,'    ',new ) ;
            end
            fc = fckeep ;
            fca = fck2 ;
            if niter == 0
              fc = fcmin ;
            end
            qgenj = 1 ;
          else
%           ------------------------------------------------------
%           4 Preparations to start the following iteration step
%           ------------------------------------------------------
%           4.1 Print values
            if mprmon >= 3  && ~ qordi
              if mprtim ~= 0
                monon(5) ;
              end
              n1prv2(dlevfn,sqrt(sumx/n),fc,niter+1, mprmon,lumon,qmixio,'*') ;
              if mprtim ~= 0
                monoff(5) ;
              end
            end
%           Print the natural level of the current iterate and return
%           it in one-step mode
            sumxs = sumx ;
            sumx = sumxa ;
            if mprsol >= 2 && niter ~= 0
              if mprtim ~= 0
                monon(5) ;
              end
              n1sout(n,xa,2,iopt,wk,mprsol,lusol) ;
              if mprtim ~= 0
                monoff(5) ;
              end
            elseif mprsol >= 1 && niter == 0
              if mprtim ~= 0
                monon(5) ;
              end
              n1sout(n,xa,1,iopt,wk,mprsol,lusol) ;
              if mprtim ~= 0
                monoff(5) ;
              end  
            end
            niter = niter+1 ;
            wk.niter = niter ;
            xiter(1:n,niter+1) = x ;
            dlevf = dlevfn ;
%     Exit Repeat If ...
            if niter >= nitmax
              ierr = 2 ;
              break ;
            end
            fckeep = fc ;
%           ------------------------------------------------------
%           4.2 Return, if in one-step mode
%
% Exit function If ...
            if mode == 1
              qsucc = 1 ;
              return ;
            end
          end
        end
%       End Repeat
%       End of main iteration loop
%       ==========================
%       ----------------------------------------------------------
%       9 Exits
%       ----------------------------------------------------------
%       9.1 Solution exit
        aprec = -1.0 ;
% 
        if ierr == 0  ||  ierr == 4
          if nonlin ~= 1
            if  ~ qordi
              if  ierr == 0 
                aprec = sqrt(sumx/n) ;
                x = x+dxq ;
                xiter(1:n,niter+2) = x ;
              else 
                aprec = sqrt(sumxa/n) ;
                if alphaa > zero  &&  iormon == 3
                  x = x+dx ;
                  xiter(1:n,niter+2) = x ;
                end
              end
%             Print final monitor output
              if mprmon >= 2
                if ierr == 0
                  if mprtim ~= 0
                    monon(5) ;
                  end
                  n1prv2(dlevfn,sqrt(sumx/n),fc,niter+1,mprmon,lumon,qmixio,'*') ;
                  if mprtim ~= 0
                    monoff(5) ;
                  end
                elseif iormon == 3
                  if mprtim ~= 0
                    monon(5) ;
                  end
                  n1prv1(dlevfn,sqrt(sumxa/n),fc,niter,new,mprmon,lumon,qmixio) ;
                  if mprtim ~= 0
                    monoff(5) ;
                  end
                end
              end
              if   iormon >= 2 
                if  iconv <= 1  &&  alphae  ~=  zero  &&  alphak  ~=  zero
                  ierr = 5 ;
                end
              end
            else
%           if qordi
              aprec = sqrt(sumxa/n) ;
            end
            if mprmon >= 1
              if qordi  ||  ierr == 4
                nout = niter ;
              else
                nout =  niter+1 ;
              end 
              fprintf(lumon,'\n\n\n %s%3i%s\n\n %s%10.3e\n', ...
                            'Solution of nonlinear system of equations obtained within ', ...
                             nout,' iteration steps', ...
                             'Achieved relative accuracy',aprec) ;
            end
          else
            if MPRMON >= 1
              fprintf(lumon,'\n\n\n%s\n\n%s\n', ...
                ' Solution of linear system of equations obtained by NLEQ1', ...
                ' No estimate available for the achieved relative accuracy')
            end
          end
        end
%       ----------------------------------------------------------
%       9.2 Fail exit messages
%       ----------------------------------------------------------
%       9.2.1 Termination, since jacobian matrix became singular
        if ierr == 1 && mprerr >= 1
          fprintf(luerr,'\n%s\n',' Iteration terminates due to singular jacobian matrix') ;
        end
%       ----------------------------------------------------------
%       9.2.2 Termination after more than NITMAX iterations
        if ierr == 2 && mprerr >= 1
          fprintf(luerr,'\n%s%3i%s\n',' Iteration terminates after nitmax ',nitmax, ...
                                      '  iteration steps' ) ;
        end
%       ----------------------------------------------------------
%       9.2.3 Damping factor FC became too small
        if ierr == 3 && mprerr >= 1
          fprintf(luerr,'\n%s%10.3e\n',' Damping factor has become too small: lambda =',fc) ;
        end
%WEI
%       ----------------------------------------------------------
%       9.2.4.1 Superlinear convergence slowed down
        if ierr == 4 && mprerr >= 1
          if iconv == 2
            ctyp = 'superlinear' ;
          end
          if iconv == 3
            ctyp = 'quadratic' ;
          end
          if qmstop
            fprintf(luerr,'\n%s%s%s\n%s\n', ...
                 ' Warning: Monotonicity test failed after ',ctyp, ...
                 ' convergence was already checked;', ...
                 ' rtol requirement may be too stringent') ;
          else
            fprintf(luerr,'\n%s%s%s\n%s\n', ...
                 ' Warning: ',ctyp,' Convergence slowed down;', ...
                 ' rtol requirement may be too stringent') ;
          end
        end
%       ----------------------------------------------------------
%       9.2.4.2 Convergence criterion satisfied before superlinear
%               convergence has been established
        if ierr == 5 && dlevfn == zero
          ierr = 0 ;
        end
        if ierr == 5 && mprerr >= 1
          fprintf(luerr,'\n%s\n%s\n%s\n', ...
                  ' Warning: No quadratic or superlinear convergence established yet', ...
                  '          your solution may perhaps may be less accurate',          ...
                  '          as indicated by the standard error estimate') ;
        end
%       ----------------------------------------------------------
%       9.2.5 Error exit due to dimension mismatch of start vector
%             and problem function output
        if ierr == 22 && mprerr >= 1
          fprintf(luerr,'Error: Dimensions of startvector and problem function output differ:\n') ;
          fprintf(luerr,'       length(x0) = %5i,   length(fcn(x0)) = %5i\n',n,length(f)) ;
        end
%       ----------------------------------------------------------
%       9.2.6 Error exit due to linear solver routine N1FACT
        if ierr == 80 && mprerr >= 1
          fprintf(luerr,'%6s%5i%s\n',' Error ',ifail, ...
                        ' signalled by linear solver N1FACT') ;
        end
%       ----------------------------------------------------------
%       9.2.7 Error exit due to linear solver routine N1SOLV
        if ierr == 81 && mprerr >= 1
          fprintf(luerr,'%6s%5i%s\n',' Error ',ifail, ...
                        ' signalled by linear solver N1SOLV') ;
        end
%       ----------------------------------------------------------
%       9.2.8 Error exit due to fail of user function FCN
        if ierr == 82 && mprerr >= 1
          fprintf(luerr,'%6s%5i%s\n',' Error ',ifail, ...
                        ' signalled by user function FCN') ;
        end
%       ----------------------------------------------------------
%       9.2.9 Error exit due to fail of user function Jacobian
        if ierr == 83 && mprerr >= 1
          fprintf(luerr,'%6s%5i%s\n',' Error ',ifail, ...
                        ' signalled by user function Jacobian') ;
        end
        if ierr >= 80 && ierr <= 83
          wk.ifail = ifail ;
        end
        if (ierr == 82 || ierr == 83) && niter <= 1 && mprerr >= 1
          fprintf(luerr,'%s\n', ...
                  ' Try to find a better initial guess for the solution') ;
        end
%       ----------------------------------------------------------
%       9.3 Common exit
        if mprerr >= 3 && ierr ~= 0 && ierr ~= 4 && nonlin ~= 1
          fprintf(luerr,'\n%s%10.3e\n', ...
                  '    Achieved relative accuracy ', conva ) ;
          aprec = conva ;
        end
        rtol = aprec ;
        sumx = sumxa ;
        if mprsol >= 2 && niter ~= 0
           mode=2 ;
           if qordi
             mode=3 ;
           end
           if mprtim ~= 0
             monon(5) ;
           end
           n1sout(n,xa,mode,iopt,wk,mprsol,lusol) ;
           if mprtim ~= 0
             monoff(5) ;
           end
        elseif mprsol >= 1 && niter == 0
           if mprtim ~= 0
             monon(5) ;
           end
           n1sout(n,xa,1,iopt,wk,mprsol,lusol) ;
           if mprtim ~= 0
             monoff(5) ;
           end
        end
        if  ~ qordi
          if ierr ~= 4
            niter = niter+1 ;
            wk.niter = niter ;
          end
          dlevf = dlevfn ;
          if mprsol >= 1
%           print solution or final iteration vector
            if ierr == 0
               modefi = 3 ;
            else
               modefi = 4 ;
            end
            if mprtim ~= 0
              monon(5) ;
            end
            n1sout(n,x,modefi,iopt,wk,mprsol,lusol) ;
            if mprtim ~= 0
              monoff(5) ;
            end
          end
        end
%       Return the latest internal scaling to XSCAL
        xscal = xw ;
%       End of exits
%       End of function N1INT
%
function xw = n1scal(n,x,xa,xscal,iscal,qinisc,iopt)
%     ------------------------------------------------------------
%
%*    Summary :
%    
%     S C A L E : To be used in connection with NLEQ1 .
%       Computation of the internal scaling vector XW used for the
%       Jacobian matrix, the iterate vector and its related
%       vectors - especially for the solution of the linear system
%       and the computations of norms to avoid numerical overflow.
%
%*    Input parameters
%     ================
%
%     N         Int     Number of unknowns
%     X(N)      Dble    Current iterate
%     XA(N)     Dble    Previous iterate
%     XSCAL(N)  Dble    User scaling passed from parameter XSCAL
%                       of interface routine NLEQ1
%     ISCAL     Int     Option ISCAL passed from IOPT-field
%                       (for details see description of IOPT-fields)
%     QINISC    Logical = 1 : Initial scaling
%                       = 0 : Subsequent scaling
%     IOPT      Struct  Options structure passed from NLEQ1 parameter list
%
%*    Output parameters
%     =================
%
%     XW(N)     Dble   Scaling vector computed by this routine
%                      All components must be positive. The follow-
%                      ing relationship between the original vector
%                      X and the scaled vector XSCAL holds:
%                      XSCAL(I) = X(I)/XW(I) for I=1,...N
%
%*    Functions called: D1MACH
%
%*    Machine constants used
%     ======================
%
%     SMALL
%     ------------------------------------------------------------
      small  = d1mach(6) ;
%*    Begin
      if iscal == 1
        xw = xscal ;
      else
        xw = zeros(n,1) ;
        for l1=1:n
          xw(l1)=max(xscal(l1),max((abs(x(l1))+abs(xa(l1)))*0.5,small)) ;
        end
      end
%$Test-Begin
      mprmon = iopt.mprmon ;
      if mprmon >= 6
        lumon = iopt.lumon ;
        fprintf(lumon,'\n\n%s\n%s\n', ...
                ' ++++++++++++++++++++++++++++++++++++++++++' , ...
                '      x-components   scaling-components    '        ) ;
        for l1 =1:n
          fprintf(lumon, '  %18.10e  %18.10e\n', x(l1), xw(l1)) ;
        end
        fprintf(lumon,'%s\n\n', ' ++++++++++++++++++++++++++++++++++++++++++' ) ;
      end
%$Test-End
%     End of function N1SCAL
%
function [aout,fw] = n1scrf(m,n,a)
%     ------------------------------------------------------------
%
%*    Summary :
%
%     S C R O W F : Row Scaling of a (M,N)-matrix in full or sparse
%                   storage mode
%
%*    Input parameters (* marks inout parameters)
%     ===========================================
%
%       M           Int    Number of rows of the matrix
%       N           Int    Number of columns of the matrix
%       A(M,N)      Dble   Matrix to be scaled
%
%*    Output parameters
%     =================
%
%       AOUT(M,N)   Dble   The row-scaled matrix
%       FW(M)       Dble   Row scaling factors - FW(i) contains
%                          the factor by which the i-th row of A
%                          has been multiplied
%
%     ------------------------------------------------------------
%*    Begin
      fw = zeros(n,1) ;
      if issparse(a)
          nza = nnz(a) ;
          [row,col] = find(a ~= 0) ;
          aout = sparse(row,col,zeros(nza,1),m,n,nza) ;
          for j=1:nza
              rj = row(j) ;
              aaj = abs(a(rj,col(j))) ;
              fw(rj) = max(fw(rj),aaj) ;
          end
          for k=1:m
             if fw(k) > 0.0
                 fw(k) = 1.0/fw(k) ;
             else
                 fw(k) = 1.0 ;
             end
          end
          for j=1:nza
              aout(row(j),col(j)) = a(row(j),col(j))*fw(row(j)) ;
          end
      else
          aout = zeros(size(a)) ;
          for k=1:m
            s1 = max(abs(a(k,1:n))) ;
            if s1 > 0.0
              s1=1.0/s1 ;
              fw(k)=s1 ;
              aout(k,1:n) = a(k,1:n)*s1 ;
            else
              fw(k)=1.0 ;
              aout(k,1:n) = a(k,1:n) ;
            end
          end
      end
%     End of function n1scrf
%
function [aout,fw] = n1scrb(n,lda,ml,mu,a)
%     ------------------------------------------------------------
%
%*    Summary :
%
%     S C R O W B : Row Scaling of a (N,N)-matrix in band storage
%                   mode
%
%*    Input parameters (* marks inout parameters)
%     ===========================================
%
%       N           Int    Number of rows and columns of the matrix
%       LDA         Int    Leading dimension of the matrix array
%       ML          Int    Lower bandwidth of the matrix (see IOPT)
%       MU          Int    Upper bandwidth of the matrix (see IOPT)
%       A(LDA,N)    Dble   Matrix to be scaled (see Note in routine
%                          DGBFA for details on storage technique)
%
%*    Output parameters
%     =================
%
%       AOUT(LDA,N) Dble   The row-scaled matrix
%       FW(N)       Dble   Row scaling factors - FW(i) contains
%                          the factor by which the i-th row of
%                          the matrix has been multiplied
%
%     ------------------------------------------------------------
%*    Begin
      fw = zeros(n,1) ;
      aout = zeros(size(a)) ;
      m2=ml+mu+1 ;
      for k=1:n
        s1=0.0 ;
        l2=max(1,k-ml) ;
        l3=min(n,k+mu) ;
        k1=m2+k ;
        for l1=l2:l3
          s1 = max(s1,abs(a(k1-l1,l1))) ;
        end
        if s1 > 0.0
          s1=1.0/s1 ;
          fw(k)=s1 ;
          for l1=l2:l3
            aout(k1-l1,l1)=a(k1-l1,l1)*s1 ;
          end
        else
          fw(k)=1.0 ;
        end
      end
%     End of function N1SCRB
%
function [l,u,p,ifail] = n1fact(n,lda,ml,mu,a,iopt)
%     ------------------------------------------------------------
%
%*    Summary :
%
%     F A C T : Call linear algebra subprogram for factorization of
%               a (N,N)-matrix
%
%*    Input parameters (* marks inout parameters)
%     ===========================================
%
%     N             Int    Order of the linear system
%     LDA           Int    Leading dimension of the matrix array A
%     ML            Int    Lower bandwidth of the matrix (only for
%                          banded systems)
%     MU            Int    Upper bandwidth of the matrix (only for
%                          banded systems)
%     A(LDA,N)      Dble   Matrix storage. See main routine NLEQ1,
%                          note 4 for how to store a banded Jacobian
%     IOPT(50)      Int    Option vector passed from NLEQ1
%
%*    Output parameters
%     =================
%
%     L             Dble   The lower triangle of the LU-decomposition
%                          of the matrix A stored in full mode - or -
%                          The complete LU-decomposition of the matrix A
%                          stored in band mode
%     U             Dble   The upper triangle of the LU-decomposition
%                          of the matrix A stored in full mode - or -
%                          not used for a band matrix A.
%     P             Int    The permutation matrix from the LU-decomposition
%                          of the matrix A stored in full mode - or -
%                          the permutation vector from the LU-decomposition
%                          of the matrix stored in band mode.
%     IFAIL         Int    Error indicator returned by this routine:
%                          = 0 matrix decomposition successfull
%                          = 1 decomposition failed - 
%                              matrix is numerically singular
%
%
%*    functions called:  DGBFA
%
%     ------------------------------------------------------------
%*    Begin
	  mstor = iopt.mstor ;
	  if mstor == 0
	    try
		  [l,u,p] = lu(a) ;
		  ifail = 0 ;
		catch
		  ifail = 1 ;
		end
	  elseif mstor == 1
%       band mode: l holds the complete lu-factorization of p*a
		[l,p,ifail] = dgbfa(a,lda,n,ml,mu) ;
		u = [] ;
	  end
	  if ifail ~= 0
		ifail = 1 ;
	  end
%
function [x,ifail] = n1solv(n,lda,ml,mu,l,u,p,b,iopt)
%     ------------------------------------------------------------
%
%*    Summary :
%
%     S O L V E : Call linear algebra subprogram for solution of
%                 the linear system A*Z = B
%
%*    Input parameters:
%     =================
%
%     N,LDA,ML,MU,A,IOPT,IFAIL :
%                        See description for function N1FACT.          
%     B(N)       Dble    Right hand side of the linear system
%
%*    Output parameters:
%     ==================
%
%     X(N)       Dble    Solution of the linear system
%     IFAIL      Int     Error indicator returned by this routine.
%
%     functions called: DGBSL
%
%     ------------------------------------------------------------
%*    Begin
      mstor = iopt.mstor ;
      if mstor == 0
        x = p*b ;  x = l\x ;  x=u\x ;
      elseif mstor == 1
        x = dgbsl(l,lda,n,ml,mu,p,b,0) ;
      end
      ifail = 0 ;
%
function [dxq,conv,sumx,dlevf] = n1lvls(n,dxq,dx1,xw,f,mprmon,qdscal)
%
%     ------------------------------------------------------------
%
%*    Summary :
%
%     L E V E L S : To be used in connection with NLEQ1 .
%     provides descaled solution, error norm and level functions
%
%*    Input parameters (* marks inout parameters)
%     ===========================================
%
%       N              Int    Number of parameters to be estimated
%       DXQ(N)         Dble   Array containing the descaled Newton
%                             correction, if QDSCAL = 0
%       DX1(N)         Dble   array containing the scaled Newton
%                             correction
%       XW(N)          Dble   Array containing the scaling values
%       F(N)           Dble   Array containing the residuum
%       MPRMON         Int    Print information parameter (see
%                             driver routine NLEQ1 )
%       QDSCAL         Logic  1, if descaling of DX1 required,
%                             else 0
%
%*    Output parameters
%     =================
%
%       DXQ(N)         Dble   Array containing the descaled Newton
%                             correction
%       CONV           Dble   Scaled maximum norm of the Newton
%                             correction
%       SUMX           Dble   Scaled natural level function value
%       DLEVF          Dble   Standard level function value
%
%     ------------------------------------------------------------
%*    Begin
      if qdscal
%       ------------------------------------------------------------
%       1.2 Descaling of solution DX1 ( stored to DXQ )
        dxq = dx1 .* xw ;
      end
%     ------------------------------------------------------------
%     2 Evaluation of scaled natural level function SUMX and
%       scaled maximum error norm CONV
      conv = max(abs(dx1)) ;
      sumx = sum(dx1.^2) ;
%     ------------------------------------------------------------
%     3 Evaluation of (scaled) standard level function DLEVF
      dlevf = sqrt(sum(f.^2)/n) ;
%     End of function N1LVLS
%
function [a,nfcnnew,ifail] = n1jac (fcn, par, n, lda, x, fx, yscal, ajdel, ajmin, nfcn)
%  ---------------------------------------------------------------------
%
%* Title
%
%  Evaluation of a dense Jacobian matrix using finite difference
%  approximation adapted for use in nonlinear systems solver NLEQ1
%
%* Environment       Matlab 6.0
%                    Sparc, Solaris
%* Latest Revision   July 2001
%
%
%* Parameter list description
%  --------------------------
%
%* External functions (to be supplied by the user)
%  -------------------------------------------------
%
%     [FOUT,FAIL] =FCN(X,FLAG,PAR) String  Name of problem-function
%       For a description refer to the driver function nleq1.
%
%* Input parameters (* marks inout parameters)
%  ----------------
%
%  FCN        String  Name of the problem-function
%  PAR        AnyType User-parameter, to be passed to the problem-function 
%  N          Int     Number of rows and columns of the Jacobian
%  LDA        Int     Leading Dimension of array A
%  X(N)       Dble    Array containing the current scaled
%                     iterate
%  FX(N)      Dble    Array containing FCN(X)
%  YSCAL(N)   Dble    Array containing the scaling factors
%  AJDEL      Dble    Perturbation of component k: abs(Y(k))*AJDEL
%  AJMIN      Dble    Minimum perturbation is AJMIN*AJDEL
%  NFCN       Int  *  FCN - evaluation count
%
%* Output parameters (* marks inout parameters)
%  -----------------
%
%  A(LDA,N)   Dble    Array to contain the approximated
%                     Jacobian matrix ( dF(i)/dx(j)in A(i,j))
%  NFCNNEW    Int  *  FCN - evaluation count adjusted
%  IFAIL      Int     Return code non-zero if Jacobian could not
%                     be computed.
%  ---------------------------------------------------------------------
%
%* Begin
%
      ifail = 0 ;
      nfcnnew = nfcn ;
      for k = 1:n
         w = x(k) ;
         
         su = sign(x(k)) ;
         if su == 0
           su = 1 ;
         end
         u = max(max(abs(x(k)),ajmin),yscal(k))*ajdel*su ;
         x(k) = w + u ;
%
         [fu,ifail] = feval(fcn,x,'',par) ;
         nfcnnew = nfcnnew + 1 ;
         if ifail  ~=  0
           break ;
         end
%
         x(k) = w ;
         a(1:n,k) = (fu-fx) / u ;
      end
%* End of N1JAC

function [a,nfcnnew,ifail] = n1jacb (fcn, par, n, lda, ml, x, fx, yscal, ajdel, ajmin, nfcn)
%  ---------------------------------------------------------------------
%
%* Title
%
%  Evaluation of a banded Jacobian matrix using finite difference
%  approximation adapted for use in nonlinear systems solver NLEQ1
%
%* Environment       Matlab 6.0
%                    Sparc, Solaris
%* Latest Revision   July 2001
%
%
%* Parameter list description
%  --------------------------
%
%* External functions (to be supplied by the user)
%  -------------------------------------------------
%
%     [FOUT,FAIL] =FCN(X,FLAG,PAR) String  Name of problem-function
%       For a description refer to the driver function nleq1.
%
%* Input parameters (* marks inout parameters)
%  ----------------
%
%  FCN        String  Name of the problem-function
%  PAR        AnyType User-parameter, to be passed to the problem-function 
%  N          Int     Number of rows and columns of the Jacobian
%  LDA        Int     Leading Dimension of array A
%  ML         Int     Lower bandwidth of the Jacobian matrix
%  X(N)       Dble    Array containing the current scaled
%                     iterate
%  FX(N)      Dble    Array containing FCN(X)
%  YSCAL(N)   Dble    Array containing the scaling factors
%  AJDEL      Dble    Perturbation of component k: abs(Y(k))*AJDEL
%  AJMIN      Dble    Minimum perturbation is AJMIN*AJDEL
%  NFCN       Int  *  FCN - evaluation count
%
%* Output parameters (* marks inout parameters)
%  -----------------
%
%  A(LDA,N)   Dble    Array to contain the approximated
%                     Jacobian matrix ( dF(i)/dx(j)in A(i,j))
%  NFCNNEW    Int  *  FCN - evaluation count adjusted
%  IFAIL      Int     Return code non-zero if Jacobian could not
%                     be computed.
%
%  ---------------------------------------------------------------------
%* Begin
%
      nfcnnew = nfcn ;
      ifail = 0 ;
      mu = lda - 2*ml - 1 ;
      ldab = ml + mu + 1 ;
      a = zeros(ldab,n) ;
      w = zeros(n,1) ;
      u = zeros(n,1) ;
%
      for jj = 1:ldab
         for k = jj:ldab:n
            w(k) = x(k) ;
            su = sign(x(k)) ;
            if su == 0
              su = 1 ;
            end
            u(k) = max(max(abs(x(k)),ajmin),yscal(k))*ajdel*su ;
            x(k) = w(k) + u(k) ;
         end
%
         [fu,ifail] = feval(fcn,x,'',par) ;
         nfcnnew = nfcnnew + 1 ;
         if ifail  ~=  0
           break ;
         end 
%
         for k = jj:ldab:n
            x(k) = w(k) ;
            i1 = max (1, k - mu) ;
            i2 = min (n, k + ml) ;
            mh = mu + 1 - k ;
            a(mh+i1:mh+i2,k) = (fu(i1:i2) - fx(i1:i2)) / u(k) ;
         end
%
      end
%* End of N1JACB
%

function [a,eta,nfcnnew,ifail] = n1jcf(fcn,par,n,lda,x,fx,yscal,eta,etamin,etamax,etadif,conv,nfcn)
%
%  ---------------------------------------------------------------------
%
%* Title
%
%  Approximation of dense Jacobian matrix for nonlinear systems solver
%  NLEQ1 with feed-back control of discretization and rounding errors
%
%* Environment       Matlab 6.0
%                    Sparc, Solaris
%* Latest Revision   July 2001
%
%
%* Parameter list description
%  --------------------------
%
%* External functions (to be supplied by the user)
%  -------------------------------------------------
%
%     [FOUT,FAIL] =FCN(X,FLAG,PAR) String  Name of problem-function
%       For a description refer to the driver function nleq1.
%
%* Input parameters (* marks inout parameters)
%  ----------------
%
%  FCN        String  Name of the problem-function
%  PAR        AnyType User-parameter, to be passed to the problem-function 
%  N          Int     Number of rows and columns of the Jacobian
%  LDA        Int     Leading dimension of A (LDA  >=  N)
%  X(N)       Dble    Array containing the current scaled
%                     iterate
%  FX(N)      Dble    Array containing FCN(X)
%  YSCAL(N)   Dble    Array containing the scaling factors
%  ETA(N)     Dble *  Array containing the scaled denominator
%                     differences
%  ETAMIN     Dble    Minimum allowed scaled denominator
%  ETAMAX     Dble    Maximum allowed scaled denominator
%  ETADIF     Dble    sqrt (1.1*EPMACH)
%                     EPMACH = machine precision
%  CONV       Dble    Maximum norm of last (unrelaxed) Newton correction
%  NFCN       Int  *  FCN - evaluation count
%
%* Output parameters (* marks inout parameters)
%  -----------------
%
%  A(LDA,N)   Dble    Array to contain the approximated
%                     Jacobian matrix ( dF(i)/dx(j)in A(i,j))
%  ETA(N)     Dble *  Scaled denominator differences adjusted
%  NFCNNEW    Int  *  FCN - evaluation count adjusted
%  IFAIL      Int     Return code non-zero if Jacobian could not
%                     be computed.
%
%* Constants
%  ---------
      small2 = 0.1 ;
%
%  ---------------------------------------------------------------------
%* Begin
%
      ifail = 0 ;
      nfcnnew = nfcn ;
      a = zeros(n,n) ;
      for k = 1:n
         is = 0 ;
         qfine = 0 ;
         qexit = 0 ;
         while ~ qfine
            w = x(k) ;
            su = sign(x(k)) ;
            if su == 0
              su = 1 ;
            end
            u = eta(k)*yscal(k)*su ;
            x(k) = w + u ;
            [fu,ifail] = feval(fcn,x,'',par) ;
            nfcnnew = nfcnnew + 1 ;
%           exit, if ...
            if ifail  ~=  0
              qexit = 1 ;
              break ;
            end
            x(k) = w ;
            sumd = 0.0 ;
            for i = 1:n
               hg = max (abs(fx(i)), abs(fu(i))) ;
               fhi = fu(i) - fx(i) ;
               if hg  ~=  0.0
                 sumd = sumd + (fhi/hg)^2 ;
               end
               a(i,k) = fhi / u ;
            end
            sumd = sqrt (sumd/n) ;
            qfine = 1 ;
            if sumd  ~=  0.0  &&  is  ==  0
               eta(k) = min( etamax , max(etamin , sqrt(etadif/sumd)*eta(k)) ) ;
               is = 1 ;
               qfine = conv  <  small2  ||  sumd  >=  etamin ;
            end
         end
         if qexit
           break ;
         end
      end
%* End of function N1JCF

function [a,eta,nfcnnew,ifail] = n1jcfb(fcn,par,n,lda,ml,x,fx,yscal,eta,etamin,etamax,etadif, ...
                                           conv,nfcn)
%     
%     ---------------------------------------------------------------------
%     
%* Title
%  
%  Approximation of banded Jacobian matrix for nonlinear systems solver
%  NLEQ1 with feed-back control of discretization and rounding errors
%  
%* Environment       Matlab 6.0
%                    Sparc, Solaris
%* Latest Revision   July 2001
%  
%  
%* Parameter list description
%  --------------------------
%  
%* External functions (to be supplied by the user)
%  -------------------------------------------------
%  
%     [FOUT,FAIL] =FCN(X,FLAG,PAR) String  Name of problem-function
%       For a description refer to the driver function nleq1.
%
%* Input parameters (* marks inout parameters)
%  ----------------
%
%  FCN        String  Name of the problem-function
%  PAR        AnyType User-parameter, to be passed to the problem-function 
%  N          Int     Number of rows and columns of the Jacobian
%  LDA        Int     Leading dimension of A (LDA  >=  ML+MU+1)
%  ML         Int     Lower bandwidth of Jacobian matrix
%  X(N)       Dble    Array containing the current scaled
%                     iterate
%  FX(N)      Dble    Array containing FCN(X)
%  YSCAL(N)   Dble    Array containing the scaling factors
%  ETA(N)     Dble *  Array containing the scaled denominator
%                     differences
%  ETAMIN     Dble    Minimum allowed scaled denominator
%  ETAMAX     Dble    Maximum allowed scaled denominator
%  ETADIF     Dble    sqrt (1.1*EPMACH)
%                     EPMACH = machine precision
%  CONV       Dble    Maximum norm of last (unrelaxed) Newton correction
%  NFCN       Int  *  FCN - evaluation count
%  
%* Output parameters (* marks inout parameters)
%  -----------------
%  
%  A(LDA,N)   Dble    Array to contain the approximated
%                     Jacobian matrix ( dF(i)/dx(j)in A(i,j))
%  ETA(N)     Dble *  Scaled denominator differences adjusted
%  NFCNNEW    Int  *  FCN - evaluation count adjusted
%  IFAIL      Int     Return code non-zero if Jacobian could not
%                     be computed.
%  
%* Constants
%  ---------
      small2 = 0.1 ;
%  
%  ---------------------------------------------------------------------
%* Begin
%  
      nfcnnew = nfcn ;
      mu = lda - 2*ml - 1 ;
      ldab = ml + mu + 1 ;
      a = zeros(ldab,n) ;
      w = zeros(n,1) ;
      u = zeros(n,1) ;
%  
      for jj = 1:ldab
         is = 0 ;
         qfine = 0 ;
         qexit = 0 ;
         while ~ qfine
            for  k = jj:ldab:n
               w(k) = x(k) ;
               su = sign(x(k)) ;
               if su == 0
                 su = 1 ;
               end
               u(k) = eta(k)*yscal(k)*su ;
               x(k) = w(k) + u(k) ;
            end
%     
            [fu,ifail] = feval(fcn,x,'',par) ;
            nfcnnew = nfcnnew + 1 ;
%           exit, if ...
            if ifail  ~=  0
              qexit = 1 ;
              break ;
            end
%     
            for k = jj:ldab:n
               x(k) = w(k) ;
               sumd = 0.0 ;
               i1 = max(1, k - mu) ;
               i2 = min(n, k + ml) ;
               mh = mu + 1 - k ; 
               for i = i1:i2
                  hg = max( abs(fx(i)),abs(fu(i)) ) ;
                  fhi = fu(i) - fx(i) ;
                  if hg  ~=  0.0
                    sumd = sumd + (fhi/hg)^2 ;
                  end
                  a(mh+i,k) = fhi / u(k) ;
               end
               sumd = sqrt (sumd/n) ;
               qfine = 1 ;
               if sumd  ~=  0.0  &&  is  ==  0
                  eta(k) = min ( etamax , max(etamin , sqrt (etadif / sumd)*eta(k)) ) ;
                  is = 1 ;
                  qfine = conv  <  small2  ||  sumd  >=  etamin ;
               end
            end
         end
         if qexit
           break ;
         end
      end
%* end of function n1jcfb

function dummy = n1prv1(dlevf,dlevx,fc,niter,new,mprmon,lumon,qmixio)
%     ------------------------------------------------------------
%
%*    Summary :
%
%     N 1 P R V 1 : Printing of intermediate values (Type 1 routine)
%
%*    Parameters
%     ==========
%
%     DLEVF, DLEVX   See descr. of internal double variables of N1INT
%     FC,NITER,NEW,MPRMON,LUMON
%                  See parameter description of function N1INT
%     QMIXIO Logical  = 1 , if LUMON == LUSOL
%                     = 0 , if LUMON ~= LUSOL
%
%     ------------------------------------------------------------
%*    End Prologue
%     Print Standard - and natural level
      if qmixio
         fprintf(lumon,'%s\n', ...
                 '  ******************************************************************') ;
         if mprmon >= 3
            fprintf(lumon,'%s\n', ...
                  '        It       Normf           Normx                     New'        ) ;
         end
         if mprmon == 2
            fprintf(lumon,'%s\n', ...
                  '        It       Normf           Normx         Damp.Fct.   New'        ) ;
         end
      end
      if mprmon >= 3 || niter == 0 
        fprintf(lumon,'      %4i     %10.3e      %10.3e                 %2i\n', ...
                niter,dlevf,dlevx,new) ;
      end
      if mprmon == 2 && niter ~= 0
        fprintf(lumon,'      %4i     %10.3e      %10.3e      %7.5f    %2i\n', ...
                niter,dlevf,dlevx,fc,new ) ;
      end
      if qmixio
         fprintf(lumon,'%s\n', ...
                 '  ******************************************************************') ;
      end
      dummy = 1 ;
%     End of function N1PRV1

function dummy = n1prv2(dlevf,dlevx,fc,niter,mprmon,lumon,qmixio,cmark)
%     ------------------------------------------------------------
%
%*    Summary :
%
%     N 1 P R V 2 : Printing of intermediate values (Type 2 routine)
%
%*    Parameters
%     ==========
%
%     DLEVF, DLEVX   See descr. of internal double variables of N1INT
%     FC,NITER,MPRMON,LUMON
%                  See parameter description of function N1INT
%     QMIXIO Logical  = 1 , if LUMON == LUSOL
%                     = 0 , if LUMON ~= LUSOL
%     CMARK Char*1    Marker character to be printed before DLEVX
%
%     ------------------------------------------------------------
%*    End Prologue
%     Print Standard - and natural level, and damping
%     factor
      if qmixio
         fprintf(lumon,'%s\n', ...
                 '  ******************************************************************') ;
         fprintf(lumon,'%s\n','        It       Normf           Normx         Damp.Fct.') ;
      end
      fprintf(lumon,'      %4i     %10.3e    %1s %10.3e      %7.5f\n', ...
              niter,dlevf,cmark,dlevx,fc) ;
      if qmixio
         fprintf(lumon,'%s\n', ...
                 '  ******************************************************************') ;
      end
      dummy = 1 ;
%     End of function N1PRV2

function dummy = n1sout(n,x,mode,iopt,wk,mprint,luout)
%     ------------------------------------------------------------
%
%*    Summary :
%
%     S O L O U T : Printing of iterate (user customizable routine)
%
%*    Input parameters
%     ================
%
%     N         Int Number of equations/unknowns
%     X(N)   Double iterate vector
%     MODE          =1 This routine is called before the first
%                      Newton iteration step
%                   =2 This routine is called with an intermedi-
%                      ate iterate X(N)
%                   =3 This is the last call with the solution
%                      vector X(N)
%                   =4 This is the last call with the final, but
%                      not solution vector X(N)
%     IOPT(50)  Int The option array as passed to the driver
%                   routine(elements 46 to 50 may be used
%                   for user options)
%     MPRINT    Int Solution print level 
%                   (see description of IOPT-field MPRINT)
%     LUOUT     Int the solution print unit 
%                   (see description of see IOPT-field LUSOL)
%
%
%*    Workspace parameters
%     ====================
%
%     WK    see description in driver routine
%
%     ------------------------------------------------------------
%*    Begin
      qnorm = 1 ;
      if qnorm
         if mode == 1
           fprintf(luout,'%s\n%s%5i\n\n%s\n','  Start data:','  N =',n, ...
                  '  Format: iteration-number, (x(i),i=1,...N) , Normf , Normx ') ;
           fprintf(luout,'%s\n','  Initial data:') ;
         elseif mode == 3
           fprintf(luout,'%s\n','  Solution data:') ;
         elseif mode == 4
           fprintf(luout,'%s\n','  Final data:') ;
         end
         fprintf(luout,' %5i\n',wk.niter) ;
         l2 = 0 ;
         for l1=1:n
           fprintf(luout,'%18.10e ',x(l1)) ;
           l2 =l2+1 ;
           if l2 == 3
             fprintf(luout,'%1s\n',' ') ;
             l2 = 0 ;
           end
         end
         fprintf(luout,'%18.10e %18.10e \n', wk.dlevf, sqrt(wk.sumx/n) ) ;
         if mode == 1 && mprint >= 2
           fprintf(luout,'%s\n','  Intermediate data:') ;
         elseif mode >= 3
           fprintf(luout,'%s\n','  End data:') ;
         end
      end
      dummy = 1 ;
%     End of function N1SOUT

function norm = wnorm(n,z,xw)
%     ------------------------------------------------------------
%
%*    Summary :
%
%     E N O R M : Return the norm to be used in exit (termination)
%                 criteria
%
%*    Input parameters
%     ================
%
%     N         Int Number of equations/unknowns
%     Z(N)     Dble  The vector, of which the norm is to be computed
%     XW(N)    Dble  The scaling values of Z(N)
%
%*    Output
%     ======
%
%     WNORM(N,Z,XW)  Dble  The mean square root norm of Z(N) subject
%                          to the scaling values in XW(N):
%                          = Sqrt( Sum(1,...N)((Z(I)/XW(I))**2) / N )
%
%     ------------------------------------------------------------
%*    Begin
      norm = sqrt( sum( (z./xw).^2 ) / n ) ;
%     End of function WNORM
%*    End package
%
%*    Group  Linear Solver functions from LINPACK (incl. BLAS)
%
function [abd,ipvt,info] = dgbfa(abd,lda,n,ml,mu)
%
%     dgbfa factors a double precision band matrix by elimination.
%
%     dgbfa is usually called by dgbco, but it can be called
%     directly with a saving in time if  rcond  is not needed.
%
%     on entry
%
%        abd     double precision(lda, n)
%                contains the matrix in band storage.  the columns
%                of the matrix are stored in the columns of  abd  and
%                the diagonals of the matrix are stored in rows
%                ml+1 through 2*ml+mu+1 of  abd .
%                see the comments below for details.
%
%        lda     integer
%                the leading dimension of the array  abd .
%                lda must be  >=  2*ml + mu + 1 .
%
%        n       integer
%                the order of the original matrix.
%
%        ml      integer
%                number of diagonals below the main diagonal.
%                0  <=  ml  <  n .
%
%        mu      integer
%                number of diagonals above the main diagonal.
%                0  <=  mu  <  n .
%                more efficient if  ml  <=  mu .
%     on return
%
%        abd     an upper triangular matrix in band storage and
%                the multipliers which were used to obtain it.
%                the factorization can be written  a = l*u  where
%                l  is a product of permutation and unit lower
%                triangular matrices and  u  is upper triangular.
%
%        ipvt    integer(n)
%                an integer vector of pivot indices.
%
%        info    integer
%                = 0  normal value.
%                = k  if  u(k,k)  ==  0.0 .  this is not an error
%                     condition for this function, but it does
%                     indicate that dgbsl will divide by zero if
%                     called.  use  rcond  in dgbco for a reliable
%                     indication of singularity.
%
%     band storage
%
%           if  a  is a band matrix, the following program segment
%           will set up the input.
%
%                   ml = (band width below the diagonal) ;
%                   mu = (band width above the diagonal) ;
%                   m = ml + mu + 1 ;
%                   for j = 1:n
%                      i1 = max(1, j-mu) ;
%                      i2 = min(n, j+ml) ;
%                      for i = i1:i2
%                         k = i - j + m ;
%                         abd(k,j) = a(i,j) ;
%                      end
%                   end
%
%           this uses rows  ml+1  through  2*ml+mu+1  of  abd .
%           in addition, the first  ml  rows in  abd  are used for
%           elements generated during the triangularization.
%           the total number of rows needed in  abd  is  2*ml+mu+1 .
%           the  ml+mu by ml+mu  upper left triangle and the
%           ml by ml  lower right triangle are not referenced.
%
%     linpack. this version dated 08/14/78 .
%     cleve moler, university of new mexico, argonne national lab.
%
%     functions
%
%     blas daxpy,dscal,idamax
      m = ml + mu + 1 ;
      info = 0 ;
      ipvt = zeros(n,1);
%
%     zero initial fill-in columns
%
      j0 = mu + 2 ;
      j1 = min(n,m) - 1 ;
      if j1  >=  j0
		for jz = j0:j1
		   i0 = m + 1 - jz ;
		   abd(i0:ml,jz) = zeros(ml-i0+1,1) ;
	    end
      end
      jz = j1 ;
      ju = 0 ;
%
%     gaussian elimination with partial pivoting
%
      nm1 = n - 1 ;
      if nm1  >= 1
		for k = 1:nm1
		   kp1 = k + 1 ;
  %
  %        zero next fill-in column
  %
		   jz = jz + 1 ;
		   if jz  <=  n && ml >= 1
		      abd(1:ml,jz) = zeros(ml,1) ;
	       end
  %
  %        find l = pivot index
  %
		   lm = min(ml,n-k) ;
		   l = idamax(lm+1,abd(m:m+lm,k),1) + m - 1 ;
		   ipvt(k) = l + k - m ;
  %
  %        zero pivot implies this column already triangularized
  %
		   if abd(l,k)  ~=  0.0
  %
  %           interchange if necessary
  %
			  if l  ~=  m
				 t = abd(l,k) ;
				 abd(l,k) = abd(m,k) ;
				 abd(m,k) = t ;
	          end
  %
  %           compute multipliers
  %
			  t = -1.0/abd(m,k) ;
			  abd(m+1:m+lm,k) = dscal(lm,t,abd(m+1:m+lm,k),1) ;
  %
  %           row elimination with column indexing
  %
			  ju = min(max(ju,mu+ipvt(k)),n) ;
			  mm = m ;
			  if ju  >=  kp1
                for j = kp1:ju
                   l = l - 1 ;
                   mm = mm - 1 ;
                   t = abd(l,j) ;
                   if l  ~=  mm
                      abd(l,j) = abd(mm,j) ;
                      abd(mm,j) = t ;
                   end
                   abd(mm+1:mm+lm,j) = daxpy(lm,t,abd(m+1:m+lm,k),1,abd(mm+1:mm+lm,j),1) ;
                end
	          end
	       else
			  info = k ;
		   end
	    end
      end
      ipvt(n) = n ;
      if abd(m,n)  ==  0.0 
        info = n ;
      end

function x = dgbsl(abd,lda,n,ml,mu,ipvt,b,job)
%
%     dgbsl solves the double precision band system
%     a * x = b  or  trans(a) * x = b
%     using the factors computed by dgbco or dgbfa.
%
%     on entry
%
%        abd     double precision(lda, n)
%                the output from dgbco or dgbfa.
%
%        lda     integer
%                the leading dimension of the array  abd .
%
%        n       integer
%                the order of the original matrix.
%
%        ml      integer
%                number of diagonals below the main diagonal.
%
%        mu      integer
%                number of diagonals above the main diagonal.
%
%        ipvt    integer(n)
%                the pivot vector from dgbco or dgbfa.
%
%        b       double precision(n)
%                the right hand side vector.
%
%        job     integer
%                = 0         to solve  a*x = b ,
%                = nonzero   to solve  trans(a)*x = b , where
%                            trans(a)  is the transpose.
%
%     on return
%
%        b       the solution vector  x .
%
%     error condition
%
%        a division by zero will occur if the input factor contains a
%        zero on the diagonal.  technically this indicates singularity
%        but it is often caused by improper arguments or improper
%        setting of lda .  it will not occur if the functions are
%        called correctly and if dgbco has set rcond  >  0.0
%        or dgbfa has set info  ==  0 .
%
%     to compute  inverse(a) * c  where  c  is a matrix
%     with  p  columns
%           call dgbco(abd,lda,n,ml,mu,ipvt,rcond,z)
%           if rcond is too small) go to ...
%           do 10 j = 1, p
%              call dgbsl(abd,lda,n,ml,mu,ipvt,c(1,j),0)
%        10 continue
%
%     linpack. this version dated 08/14/78 .
%     cleve moler, university of new mexico, argonne national lab.
%
%     functions
%
%     blas daxpy,ddot
%
      m = mu + ml + 1 ;
      nm1 = n - 1 ;
      if job  ==  0
%
%        job = 0 , solve  a * x = b
%        first solve l*y = b
%
         if ml  ~=  0 && nm1 >= 1
            for k = 1:nm1
               lm = min(ml,n-k) ;
               l = ipvt(k) ;
               t = b(l) ;
               if l  ~=  k 
                  b(l) = b(k) ;
                  b(k) = t ;
               end
               b(k+1:k+lm) = daxpy(lm,t,abd(m+1:m+lm,k),1,b(k+1:k+lm),1) ;
            end
         end
%
%        now solve  u*x = y
%
         for kb = 1:n
            k = n + 1 - kb ;
            b(k) = b(k)/abd(m,k) ;
            lm = min(k,m) - 1 ;
            la = m - lm ;
            lb = k - lm ;
            t = -b(k) ; 
            b(lb:lb+lm-1) = daxpy(lm,t,abd(la:la+lm-1,k),1,b(lb:lb+lm-1),1) ;
         end
      else
%
%        job = nonzero, solve  trans(a) * x = b
%        first solve  trans(u)*y = b
%
         for k = 1:n
            lm = min(k,m) - 1 ;
            la = m - lm ;
            lb = k - lm ;
            t = ddot(lm,abd(la:la+lm-1,k),1,b(lb:lb+lm-1),1) ;
            b(k) = (b(k) - t)/abd(m,k) ;
         end
%
%        now solve trans(l)*x = y
%
         if ml  ~=  0 && nm1 >= 1
            for kb = 1:nm1
               k = n - kb ;
               lm = min(ml,n-k) ;
               b(k) = b(k) + ddot(lm,abd(m+1:m+lm,k),1,b(k+1:k+lm),1) ;
               l = ipvt(k) ;
               if l  ~=  k
                  t = b(l) ;
                  b(l) = b(k) ;
                  b(k) = t ;
               end
            end
         end
      end
      x = b ;

function dyout = daxpy(n,da,dx,incx,dy,incy)
%
%     constant times a vector plus a vector.
%     jack dongarra, linpack, 3/11/78.
%
      if n > 0
        if incx ~= 1
          dyout = zeros(size(dy)) ;
        end
        if incx > 0 && incy > 0
           dyout(1:incy:(n-1)*incy+1) = dy(1:incy:(n-1)*incy+1) + da * dx(1:incx:(n-1)*incx+1) ;
        else
          error('daxpy: incx and incy must be positive') ;
        end ;
      else
        dyout = dy ;
      end ;

function dxout = dscal(n,da,dx,incx)
%
%     scales a vector by a constant.
%     jack dongarra, linpack, 3/11/78.
%
      if n > 0
        if incx ~= 1
          dxout = zeros(size(dx)) ;
        end
        if incx > 0
          dxout(1:incx:n*incx) = da*dx(1:incx:n*incx) ;
        else
          error('dscal: incx must be positive') ;
        end
      end

function idamx = idamax(n,dx,incx)
%
%     finds the index of element having max. absolute value.
%     jack dongarra, linpack, 3/11/78.
%
      if n < 1
        idamx = 0 ;
        return ;
      end
      if n == 1
        idamx = 1 ;
        return ;
      end
%
%        code for general increment
%
      ix = 1 ;
      idamx = 1 ;
      dmax = abs(dx(1)) ;
      ix = ix + incx ;
      for i = 2:n
         if abs(dx(ix)) > dmax
           idamx = i ;
           dmax = abs(dx(ix)) ;
         end
         ix = ix + incx ;
      end

function dot = ddot(n,dx,incx,dy,incy)
%
%     forms the dot product of two vectors.
%     jack dongarra, linpack, 3/11/78.
%
      dot = 0.0 ;
      if N > 0
        if incx > 0 && incy > 0
          dot = sum( dx(1:incx:(n-1)*incx+1) .* dy(1:incy:(n-1)*incy+1) ) ;
        else
          error('ddot: incx and incy must be positive') ;
        end
      end

%*    Group  Time monitor package
%
%*    Begin Prologue
%     ------------------------------------------------------------
%
%*  Title
%    
%     Monitor - A package for making multiple time measurements and
%               summary statistics
%
%*  Written by        U. Nowak, L. Weimann 
%*  Version           1.0
%*  Revision          July 2001
%*  Latest Change     July 2001
%*  Library           CodeLib
%*  Code              Fortran 77, Double Precision
%*  Environment       Standard Fortran 77 environment on PC's,
%*  Copyright     (c) Konrad Zuse Zentrum fuer
%                     Informationstechnik Berlin
%                     Takustrasse 7, D-14195 Berlin-Dahlem
%                     phone : + 49/30/84185-0
%                     fax   : + 49/30/84185-125
%*  Contact           Bodo Erdmann 
%                     ZIB, Numerical Software Development 
%                     phone : + 49/30/84185-185
%                     fax   : + 49/30/84185-107
%                     e-mail: erdmann@zib.de
%
%  ---------------------------------------------------------------
%
%* Licence
%    You may use or modify this code for your own non commercial
%    purposes for an unlimited time. 
%    In any case you should not deliver this code without a special 
%    permission of ZIB.
%    In case you intend to use the code commercially, we oblige you
%    to sign an according licence agreement with ZIB.
%
%* Warranty 
%    This code has been tested up to a certain level. Defects and
%    weaknesses, which may be included in the code, do not establish
%    any warranties by ZIB. ZIB does not take over any liabilities
%    which may follow from aquisition or application of this code.
%
%* Software status 
%    This code is under care of ZIB and belongs to ZIB software class 1.
%
%  ---------------------------------------------------------------
%
%*    Summary:
%
%     Monitor is a package for generating time and summary statistics
%     about the execution of multiple program parts of any program.
%     Nested measurements of program parts are possible.
%     ------------------------------------------------------------
%
%*    Usage:
%
%     The usage of Monitor is naturally divided into three phases:
%     1. the initialization and setup phase before the start of
%        the program or functions package to be measured;
%     2. the run phase of the program to be measured;
%     3. the final evaluation call.
%
%     The phase 1 must start with exactly one call of the function
%     MONINI, which passes a title string and a logical unit for
%     later statistics output and possible error messages to the
%     package. This call follows a number of calls of the function
%     MONDEF, where each call associates an identification string
%     to a positive integer number, called the measurement index
%     - up to maxtab, where maxtab is a package constant. Multiple
%     measurement indices may be used for measurements of multiple
%     program parts. The index 0 must also be associated with some
%     identification string, and corresponds to all parts of the
%     measured program from the measurement start call till the final
%     evaluation call, which are not associated with specific positive
%     measurement indices. After all necessary MONDEF calls are done,
%     the measurements are started at begin of the program to be
%     measured by a parameterless call of MONSRT.
%     In phase 2, each program part to be measured must be immediately
%     preceeded by a call of the function MONON with the associated 
%     measurement index, and must be immediately followed by a call of
%     the function MONOFF with the same measurement index. Measure-
%     ments of nested program parts are possible, and nesting is allowed
%     up to the number mnest, where mnest is a package constant.
%     Calling MONOFF without a preceeding MONON call with the same 
%     measurement index, or calling one of these functions with a
%     measurement index not previously defined by a MONDEF call causes
%     an error stop of the program. 
%     Finally at the end of the program to be measured, a 
%     call of the function MONEND closes all measurements and
%     prints the summary statistics.
%     As delivered, maxtab has a value 20 and mnest a value 10, but
%     both constants may be increased, if needed, to any possible
%     integer value, by simply changing it's values where they
%     are set, below.
%
%*    functions and their parameters:
%     =================================
%
%     MONINI(CIDENT,LUMON)  : Initialize Monitor
%       CIDENT  char*20  Identification string for the total measurement
%                        ( printed in summary )
%       LUMON   int      The logical unit for printing out the summary
%
%     MONDEF(MESIND,CIDMES) : Define one measurement index
%       MESIND  int      >=1 : measurement index for a specific part
%                        = 0 : measurement index for all remaining parts
%                              (i.e. not belonging to parts with 
%                               index >=1)
%       CIDMES  char*15  Identification string for the part associated
%                        with MESIND ( printed in summary )
%
%     MONSRT (dummy)        : Start measurements
%
%     MONON(MESIND)         : Start measurement of a specific part
%       MESIND  int      >=1 : measurement index for a specific part
%
%     MONOFF(MESIND)        : Stop measurement of a specific part
%       MESIND  int      >=1 : measurement index for a specific part
%
%     MONEND (DUMMY)         : Finish measurements and print summary
%       
%
%
%*    Example:
%       Calling sequence:
%
%        monini (' example',1) ;
%        mondef (0,'solver') ;
%        mondef (1,'user function') ;
%        mondef (2,'user matrix') ;
%        monsrt (dummy) ;
%        ...
%        program to be measured (part without specific measurement index)
%        ...
%        while ~ termination      
%          ...
%          monon (2) ;
%          ...  user matrix code ...
%          monoff(2) ;
%          ...
%          program to be measured (part without specific measurement index)
%          ...
%          monon (1) ;
%          ...  user function code ...
%          monoff(1) ;
%          ...
%          program to be measured (part without specific measurement index)
%          ...
%        end
%        ...
%        monend (dummy) ;
%     ------------------------------------------------------------
%
%     initialize monitor
%
function dummy = monini(texth,iounit)
%
      global maxtab mnest ;
      global sec count asec pc1 pc2 indxo time1 time0 maxind name;
      global sec0 name0 text moni qon ioncnt indact info;
      maxtab = 20 ; mnest = 10 ;
      info = 1 ;
      moni=iounit ;
      maxind=0 ;
      text=texth ;
      sec(1:maxtab)=zeros(maxtab,1) ;
      asec(1:maxtab)=zeros(maxtab,1) ;
      count(1:maxtab)=zeros(maxtab,1) ;
      qon(1:maxtab)=zeros(maxtab,1) ;
      indact(1:mnest)=zeros(mnest,1) ;
      sec0=0.0 ;
      ioncnt=0 ;
      dummy = 1 ;
%
%     define one monitor entry
%
function dummy = mondef(indx,nameh)
      global maxtab  mnest ;
      global sec count asec pc1 pc2 indxo time1 time0 maxind name;
      global sec0 name0 text moni qon ioncnt indact ;
      if indx >= 0  &&  indx <= maxtab
        if indx > maxind
          maxind=indx ;
        end
		if indx > 0
		  if count(indx) > 0
		    error('\n%s\n%s%4i%s\n',' error in function mondef', ...
		          '   indx = ',indx,'   already in use') ;
          end
		end
		if indx == 0
		  name0 = nameh ;
		else
		  name(indx,1:length(nameh)) = nameh ;
		end
      else
        error('\n%s\n%s%4i\n%s%4i\n',' error in function mondef', ...
              '   indx > ',maxtab,'   indx=',indx) ;
      end
      dummy = 1 ;
%
%     start monitor measurements
% 
function dummy = monsrt(dummy2)
      global time1 ;
      time1 = cputime ;
      dummy = 1 ;
%      if igraph > 0
%        gmini(maxind,name0,name)
%      end

%
%     start one measurement
%
function dummy = monon(indx)
      global maxtab  mnest ;
      global sec count asec pc1 pc2 indxo time1 time0 maxind name;
      global sec0 name0 text moni qon ioncnt indact info;
      if indx > maxind || indx <= 0
        error('\n%s\n%s%4i\n',' error in function monon', ...
                      '   indx out of range    indx=',indx) ;
      end
      if qon(indx)
        error('\n%s\n%s%4i\n',' error in function monon', ...
              '   measurement is already running for    indx=',indx) ;
      end
      asec(indx) = cputime ;
      qon(indx)=1 ;
      if ioncnt == 0
        sec0=sec0+asec(indx)-time1 ;
      else
        indxo=indact(ioncnt) ;
        sec(indxo)=sec(indxo)+asec(indx)-asec(indxo) ;
      end
      ioncnt=ioncnt+1 ;
      indact(ioncnt)=indx ;
      if info > 1
        fprintf(moni,'%7s%15s%7.5f',' enter ',name(indx,:),asec(indx)) ;
      end
      dummy = 1 ;
%
%      if igraph > 0
%        call gmon(indx,sec0)
%      end

%
%     stop one measurement
%
function  dummy = monoff(indx)
      global maxtab  mnest ;
      global sec count asec pc1 pc2 indxo time1 time0 maxind name;
      global sec0 name0 text moni qon ioncnt indact info;
      if indx > maxind || indx <= 0
        error('\n%s\n%s%4i\n',' error in function monoff', ...
                      '   indx out of range    indx=',indx) ;
      end
      if ~ qon(indx)
        error('\n%s\n%s%4i\n',' error in function monoff', ...
              '   measurement has never been activated for     indx=',indx) ;
      end
      time2 = cputime ;
      qon(indx)=0 ;
      sec(indx)=sec(indx)+time2-asec(indx) ;
      count(indx)=count(indx)+1 ;
      ioncnt=ioncnt-1 ;
      if ioncnt == 0
        time1=time2 ;
      else
        asec(indact(ioncnt))=time2 ;
      end
      if info > 1
        fprintf(moni,'%6s%15s%7.5f',' exit ',name(indx,:),asec(indx)) ;
      end
      dummy = 1 ;
%
%      if igraph > 0
%        call gmoff(indx,sec(indx))
%      end
%

%
%     terminate monitor and print statistics
%
function dummy = monend(dummy2)
      global maxtab  mnest ;
      global sec count asec pc1 pc2 indxo time1 time0 maxind name;
      global sec0 name0 text moni qon ioncnt indact ;
      time0 = cputime ;
      sec0=sec0+time0-time1 ;
%
      sum1=1.0e-10 ;
      sum1 = sum1 + sum( sec(1:maxind)) ;
      for i=1:maxind
        if count(i) > 0
          asec(i)=sec(i)/count(i) ;
        end
      end
      sum0=sum1+sec0 ;
%
      pc1(1:maxind)=100.0*sec(1:maxind)/sum0 ;
      pc2(1:maxind)=100.0*sec(1:maxind)/sum1 ;
      pc10=100.0*sec0/sum0 ;
      pc20=100.0*sec0/sum1 ;
%
      fprintf(moni,'\n\n%1s\n',' ') ;
      fprintf(moni,'%s\n', ...
              ' ###########################################################################') ;
      fprintf(moni,'%s\n%s\n', ...
           ' #                                                                         #', ...
           ' #                                                                         #') ;
      fprintf(moni,'%s%-29s%3s\n',' #   Results from time monitor program for: ', text,'  #') ;
      fprintf(moni,'%s\n%', ...
              ' #                                                                         #' ) ;
      fprintf(moni,'%16s%11.3f     %13s%11.3f                  %1s\n', ...
                   ' #   Total time:',sum0,'Sum of parts:',sum1,' #') ;
      fprintf(moni,'%s\n%', ...
              ' #                                                                         #' ) ;
      fprintf(moni, ...
              ' #     name            calls       time    av-time    %% total      %% sum   #\n') ;
%
      i0=1 ;
      fprintf(moni,' #   %-15s%8i%11.3f%11.4f%11.2f              #\n',name0,i0,sec0,sec0,pc10) ;
%
      for i=1:maxind
        fprintf(moni,' #   %-15s%8i%11.3f%11.4f%11.2f%11.2f   #\n', ...
                     char(name(i,:)),count(i),sec(i),asec(i),pc1(i),pc2(i) ) ;
      end
%
      fprintf(moni,'%s\n%', ...
              ' #                                                                         #' ) ;
      fprintf(moni,'%s\n', ...
              ' ###########################################################################') ;
      fprintf(moni,'\n\n%1s\n',' ') ;
      dummy = 1 ;
%      if igraph > 0
%         gmend ;
%      end
%
%  end package monitor

%
%*    Group  Machine dependent functions

function dmach = d1mach(i)
%
%  DOUBLE-PRECISION MACHINE CONSTANTS
%
%  D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
%
%  D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
%
%  D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING.
%
%  D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING.
%
%  D1MACH( 5) = LOG10(B)
%
%  TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT,
%  THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY
%  REMOVING THE C FROM COLUMN 1.
%  ON RARE MACHINES A STATIC STATEMENT MAY NEED TO BE ADDED.
%  (BUT PROBABLY MORE SYSTEMS PROHIBIT IT THAN REQUIRE IT.)
%
%  FOR IEEE-ARITHMETIC MACHINES (BINARY STANDARD), ONE OF THE FIRST
%  TWO SETS OF CONSTANTS BELOW SHOULD BE APPROPRIATE.
%
%  WHERE POSSIBLE, DECIMAL, OCTAL OR HEXADECIMAL CONSTANTS ARE USED
%  TO SPECIFY THE CONSTANTS EXACTLY.  SOMETIMES THIS REQUIRES USING
%  EQUIVALENT INTEGER ARRAYS.  IF YOUR COMPILER USES HALF-WORD
%  INTEGERS BY DEFAULT (SOMETIMES CALLED INTEGER*2), YOU MAY NEED TO
%  CHANGE INTEGER TO INTEGER*4 OR OTHERWISE INSTRUCT YOUR COMPILER
%  TO USE FULL-WORD INTEGERS IN THE NEXT 5 DECLARATIONS.
%
%     dmach(1) = realmin ;
%     dmach(2) = realmax ;
%     dmach(3) = eps ;
%     dmach(4) = 2.0*eps ;
%     dmach(5) = log10(2.0) ;
      if i  <  1   ||   i  >  6
        error(' d1mach - i out of bounds%10i',i) ;
      end
      if i  ==  1 
        dmach = realmin ;
      elseif i  ==  2
        dmach = realmax ;
      elseif i  ==  3
        dmach = eps ;
      elseif i  ==  4
        dmach = 2.0*eps ;
      elseif i  ==  5
        dmach = log10(2.0) ;
      elseif i  ==  6
%       dmach = sqrt(dmach(1)/dmach(3)) ;
        dmach = 4.94d-32 ;
      end

function [optval,opt] = getopt(opt,field,default)
      if isfield(opt,field)
        optval = getfield(opt,field) ;
      else
        optval = default ;
        opt = setfield(opt,field,default) ;
      end
      
function opt = iniopt(opt,field,default)
      if ~ isfield(opt,field)
        opt = setfield(opt,field,default) ;
      end
      
