      SUBROUTINE GBIT1E(N,NELT,LA,A,LINDAW,INDAWK,RHS,X,XSCAL,RTOL,IOPT,
     $                  IERR,LRWK,RWK)
C*    Begin Prologue GBIT1E
      INTEGER N, NELT, LA
      DOUBLE PRECISION A(NELT)
      INTEGER LINDAW, INDAWK(LINDAW)
      DOUBLE PRECISION RHS(N),X(N),XSCAL(N),RTOL
      INTEGER IOPT(50),IERR
      INTEGER LRWK
      DOUBLE PRECISION RWK(LRWK)
C     ------------------------------------------------------------
C
C*  Title
C
C     Good Broyden easy to use subroutine -
C     Iterative solution of a linear system with matrix given in 
C     SLAP Triad or SLAP column format. 
C
C*  Written by        U. Nowak, L. Weimann
C*  Purpose           Iterative solution of large scale systems of 
C                     linear equations
C*  Method            Secant method Good Broyden with adapted
C                     linesearch
C*  Category          D2a. - Large Systems of Linear Equations
C*  Keywords          Linear Equations; Large Systems;
C*                    Iterative Methods
C*  Version           1.1
C*  Revision          March 1991
C*  Latest Change     March 1991
C*  Library           CodeLib
C*  Code              Fortran 77, Double Precision
C*  Environment       Standard Fortran 77 environment on PC's,
C                     workstations and hosts.
C*  Copyright     (c) Konrad-Zuse-Zentrum fuer
C                     Informationstechnik Berlin (ZIB)
C                     Takustrasse 7, D-14195 Berlin-Dahlem
C                     phone : + 49/30/84185-0
C                     fax   : + 49/30/84185-125
C*  Contact           Bodo Erdmann
C                     ZIB, Division Scientific Computing, 
C                          Department Numerical Analysis and Modelling
C                     phone : + 49/30/84185-185
C                     fax   : + 49/30/84185-107
C                     e-mail: erdmann@zib.de
C
C*  References 
C
C     /1/ P. Deuflhard, R. Freund, A. Walter:
C         Fast Secant Methods for the Iterative Solution of
C         Large Nonsymmetric Linear Systems.
C         ZIB, Preprint SC 90-5 (July 1990)
C
C     /2/ U. Nowak, L. Weimann:
C         GIANT - A Software Package for the Numerical Solution
C         of Very Large Systems of Highly Nonlinear Equations. 
C         ZIB, Technical Report TR 90-11 (December 1990)
C
C  ---------------------------------------------------------------
C
C* Licence
C    You may use or modify this code for your own non commercial
C    purposes for an unlimited time. 
C    In any case you should not deliver this code without a special 
C    permission of ZIB.
C    In case you intend to use the code commercially, we oblige you
C    to sign an according licence agreement with ZIB.
C
C* Warranty 
C    This code has been tested up to a certain level. Defects and
C    weaknesses, which may be included in the code, do not establish
C    any warranties by ZIB. ZIB does not take over any liabilities
C    which may follow from aquisition or application of this code.
C
C* Software status 
C    This code is under care of ZIB and belongs to ZIB software class 1.
C
C  ---------------------------------------------------------------
C
C
C*    Parameters list description (* marks inout parameters)
C     ======================================================
C    
C*    Input parameters (* marks inout parameters)
C     ===========================================
C
C     N          Int     Dimension of linear system
C     NELT       Int     Number of nonzero elements of the matrix
C     LA         Int     Length of the array A. The required minimum
C                        depends on the selected preconditioner
C                        ( see IPRE=IOPT(35) below ).
C                        Required minimum amount is:
C                        for IPRE=0: 2*NELT
C                        for IPRE=1: NELT + INT( (NELT+N)/2)
C                        for IPRE=2: NELT + N
C                        for IPRE=3: NELT
C                        for IPRE=4: NELT + N*NBSIZE
C                        ( for NBSIZE see IOPT(36) below)
C   * A          Dble    The first NELT positions hold the real values
C                        of the matrix stored in SLAP Triad or Column 
C                        format. The locations up from A(NELT+1) are
C                        used as workspace by the preconditioner sub-
C                        routines.
C     LINDAW     Int     Length of the array INDAWK. The required
C                        minimum depends on the selected preconditioner
C                        ( see IPRE=IOPT(35) below ).
C                        Required minimum amount is: 
C                        max( 10+2*NELT , LINDPR), where  
C                        for IPRE=0: LINDPR = 2*NELT+4*N+14
C                        for IPRE=1: LINDPR = INT((3*NELT+N)/2)+2*N+13
C                        for IPRE=2: LINDPR = NELT+N+12
C                        for IPRE=3: LINDPR = NELT+N+12
C                        for IPRE=4: LINDPR = NELT+N*(NBSIZE+1)+12
C                        ( for NBSIZE see IOPT(36) below)
C   * INDAWK     Int     The first 11 positions of INDAWK are used as
C                        workspace by EASYPACK. The following elements
C                        holds the index arrays IA and JA of the SLAP
C                        format (see Note 4 for details on the formats).
C                        If INDAWK(1) .EQ. 0 on input,
C                        IA must be stored from INDAWK(12) upto
C                        INDAWK(11+NELT), and JA must be stored from
C                        INDAWK(11+NELT+1) upto INDAWK(11+2*NELT).
C                        The locations up from INDAWK(11+2*NELT+1) are
C                        used as workspace by the preconditioner sub-
C                        routines. 
C                        If INDAWK(1)=NZMAX .NE. 0 on input,
C                        IA must be stored from INDAWK(12) upto
C                        INDAWK(11+NELT), and JA must be stored from
C                        INDAWK(11+NZMAX+1) upto INDAWK(11+NZMAX+NELT).
C                        The elements representing JA are internally
C                        restored such that the same storage scheme
C                        results as requested with INDAWK(1) .EQ. 0 .
C     RHS(N)     Dble    The right hand side vector
C   * X(N)       Dble    The iteration starting vector
C   * XSCAL(N)   Dble    Scaling vector
C   * RTOL       Dble    Required (relative) precision
C   * IOPT(50)   Int     Options and integer workspace array
C                        (see below for usage)
C
C*    Output parameters
C     =================
C
C   * X(N)       Dble    The solution or final iterate vector
C   * XSCAL(N)   Dble    Scaling vector.
C                        Modified, if it containes bad input values:
C                        If (XSCAL(I).LT. SMALL) XSCAL(I) = SMALL ,
C                        If (XSCAL(I).GT. GREAT) XSCAL(I) = GREAT .
C                        For SMALL and GREAT, see section machine
C                        constants below and regard note 2.
C   * RTOL       Dble    Achieved (relative) precision (see note 1.)
C   * IOPT(50)   Int     Options and integer workspace array
C                        (see below for usage)
C     IERR       Int     Error code - a zero return signals ok.
C                        See below for nonzero error codes.
C
C*    Workspace parameters
C     ====================
C
C     LRWK           Int    Declared dimension of real workspace
C                           Required minimum: (KMAX+5)*N+2*KMAX+50
C     RWK(LRWK)      Dble   Real Workspace
C
C     Note 1.
C        The iteration stops, if the criterium for convergence is
C        RHO*sqrt(sig_mid/N) <= RTOL*norm(X_iter)
C        is satisfied, where X_iter denotes the actual iterate X and
C        sig_mid the mean sigma value 
C        sig_mid := (sig_iter-2+2*sig_iter-1+sig_iter)/4, 
C        with the raw error estimate sig_j for the iterate j. The
C        factor RHO may be modified - see RWK(21).
C        The norm used herein is defined by
C        norm(Z) := sqrt( ( Sum(i=1,N) (Z(i)/XSCAL(i))**2 )/N ).
C
C     Note 2.
C        The machine dependent values SMALL, GREAT and EPMACH are
C        gained from calls of the machine constants function ZIBCONST.
C        As delivered, this function is adapted to use constants 
C        suitable for all machines with IEEE arithmetic. If you use
C        another type of machine, you may change ZIBCONST.
C
C*   Options IOPT:
C    =============
C
C     Pos. Name   Default  Meaning
C
C       1  QSUCC    0      Indicator for the iteration modus:
C                          =0 : A new iteration will be started
C                          =1 : A previously terminated iteration will 
C                               be continued - in general without a
C                               restart
C                               (usally with a smaller tolerance RTOL
C                                prescribed as before or a larger 
C                                maximum iteration number count).
C       2..12              Reserved
C      13  MPRMON   0      Output level of print monitor
C                          = 0 : no output will be written
C                          = 1 : only a summary output will be written
C                          > 1 : reserved for future use
C                          =-j : A summary output and additionally each 
C                                j-th iterates statistics will be 
C                                written
C      14  LUMON    6      Logical unit number for print monitor
C      15..16              Reserved
C      17  MPROPT   0      Print monitor option:
C                          = 0: Standard print monitor
C                          = 1: Test print monitor for special purposes
C      18                  Reserved
C      19  MPRTIM   0      Output level for the time monitor
C                          = 0 : no time measurement and no output
C                          = 1 : time measurement will be done and
C                                summary output will be written -
C                                regard note 3.
C      20  LUTIM    6      Logical output unit for time monitor
C      21  MPRSPE   0      Output level for special information
C                          =0: No special output to be done
C                          =1: Special information will be written
C                          The special information may be used to
C                          generate a graphic display of error- and
C                          solution norms and of restart behaviour.
C      22  LUSPE    6      Logical output unit for special information
C      23..30              Reserved
C      31  KMAX     9      Maximum number of latest iterates to be saved
C                          and used by the iterative linear solver.
C                          An input <=-2 means KMAX=0 will be used. 
C                          A "-1" input means, that there no limit 
C                          applies, e.g. the core routine will be
C                          served with the maximum possible value 
C                          allowed by the total workspace amount.
C                          A zero input means KMAX=9 (or smaller
C                          downto a reasonable minimum value,
C                          if workspace is not sufficient). 
C      32  LITMAX 100      Maximum number of iterations allowed
C      33..34              Reserved
C      35  IPRE     0      Type of preconditioner, valid values are:                
C                          0: The SLAP preconditioner ILU (Incomplete LU-
C                             decomposition) will be used,
C                          1: The SLAP preconditioner lower triangle 
C                             will be used,
C                          2: The SLAP preconditioner diagonal scaling
C                             will be used,
C                          3: No preconditioner will be used,
C                          4: The ZIB supplied preconditioner block diagonal
C                             scaling will be used.
C      36  NBSIZE   2      Only, if IPRE=4: the size of the diagonal
C                          blocks. If IPRE=4 is selected and N is not a
C                          multiple of NBSIZE, an error termination
C                          occurs.
C      37..40              Reserved
C      41  NITER  (Int-Ws) Number of iteration steps totally done
C      42  K      (Int-Ws) Number of iteration steps done since
C                          latest internal restart due to restrictions
C                          implied by KMAX, TAUMIN and TAUMAX.
C      43  NMULJ  (Int-Ws) Number of (pairwise) calls of MULJAC and
C                          PRECON   
C      44  NRW    (Int-Ws) Amount of real workspace used
C      45  NRWKFR (Int-Ws) First element of RWK which is free to be used
C                          as workspace between successive calls of 
C                          GBIT1.
C      46..50              Reserved
C
C     Note 3.
C        The integrated time monitor calls the machine dependent
C        subroutine SECOND to get the actual time stamp in form
C        of a real number (Single precision). As delivered, this
C        subroutine always return 0.0 as time stamp value. Refer
C        to the compiler- or library manual of the FORTRAN compiler
C        which you actually use to find out how to get the actual
C        time stamp on your machine.
C
C     Note 4.
C       Description of the Jacobian SLAP format
C       (as extracted from the SLAP package description):
C       
C       =================== S L A P Triad format ===================
C       This routine requires that the  matrix A be   stored in  the
C       SLAP  Triad format.  In  this format only the non-zeros  are
C       stored.  They may appear in  *ANY* order.  The user supplies
C       three arrays of  length NELT, where  NELT is  the number  of
C       non-zeros in the matrix: (IA(NELT), JA(NELT), A(NELT)).  For
C       each non-zero the user puts the row and column index of that
C       matrix element  in the IA and  JA arrays.  The  value of the
C       non-zero   matrix  element is  placed  in  the corresponding
C       location of the A array.   This is  an  extremely  easy data
C       structure to generate.  On  the  other hand it   is  not too
C       efficient on vector computers for  the iterative solution of
C       linear systems.  Hence,   SLAP changes   this  input    data
C       structure to the SLAP Column format  for  the iteration (but
C       does not change it back).
C       
C       Here is an example of the  SLAP Triad   storage format for a
C       5x5 Matrix.  Recall that the entries may appear in any order.
C
C           5x5 Matrix       SLAP Triad format for 5x5 matrix on left.
C                              1  2  3  4  5  6  7  8  9 10 11
C       |11 12  0  0 15|   A: 51 12 11 33 15 53 55 22 35 44 21
C       |21 22  0  0  0|  IA:  5  1  1  3  1  5  5  2  3  4  2
C       | 0  0 33  0 35|  JA:  1  2  1  3  5  3  5  2  5  4  1
C       | 0  0  0 44  0|
C       |51  0 53  0 55|
C       
C       =================== S L A P Column format ==================
C       This routine requires  that the  matrix  A be  stored in the
C       SLAP Column format.  In this format the non-zeros are stored
C       counting down columns (except for  the diagonal entry, which
C       must appear first in each  "column")  and are stored  in the
C       double precision array A.   In other words,  for each column
C       in the matrix put the diagonal entry in  A.  Then put in the
C       other non-zero  elements going down  the column (except  the
C       diagonal) in order.   The  IA array holds the  row index for
C       each non-zero.  The JA array holds the offsets  into the IA,
C       A arrays  for  the  beginning  of each   column.   That  is,
C       IA(JA(ICOL)),  A(JA(ICOL)) points   to the beginning  of the
C       ICOL-th   column    in    IA and   A.      IA(JA(ICOL+1)-1),
C       A(JA(ICOL+1)-1) points to  the  end of the   ICOL-th column.
C       Note that we always have  JA(N+1) = NELT+1,  where N is  the
C       number of columns in  the matrix and NELT  is the number  of
C       non-zeros in the matrix.
C       
C       Here is an example of the  SLAP Column  storage format for a
C       5x5 Matrix (in the A and IA arrays '|'  denotes the end of a 
C       column):
C       
C       5x5 Matrix      SLAP Column format for 5x5 matrix on left.
C       1  2  3    4  5    6  7    8    9 10 11
C       |11 12  0  0 15|   A: 11 21 51 | 22 12 | 33 53 | 44 | 55 15 35
C       |21 22  0  0  0|  IA:  1  2  5 |  2  1 |  3  5 |  4 |  5  1  3
C       | 0  0 33  0 35|  JA:  1  4  6    8  9   12
C       | 0  0  0 44  0|
C       |51  0 53  0 55|
C
C*    Optional REAL input/output in RWK:
C     ====================================
C
C     Pos. Name          Meaning
C   
C       1..20            Reserved
C      21  RHO    IN     Security factor in error estimate (RHO.GE.1).
C                        Default: 4.0D0
C      22  TAUMIN IN     Minimum accepted stepsize factor tk.
C                        Default  1.0D-8
C      23  TAUMAX IN     Maximum accepted stepsize factor tk.
C                        Default: 1.0D2
C      24  TAUEQU IN     Parameter for "near equal" cycle check of
C                        rejected tau values. Default: 1.0D-2
C      25..30            Reserved
C      31  DIFNRM OUT    The norm2 of the latest correction of the
C                        iterate.
C      32  SOLNRM OUT    The norm2 of the final iterate.
C      33..40            Reserved
C      41  SIGMAK INTERN Norm2-square of latest raw correction DELTAK.
C      42  TAUN1  INTERN First kept tau value for tau-cycle check.
C      43  TAUN2  INTERN Second kept tau value for tau-cycle check.
C      44  SOLNRA INTERN Saves norm2 of the previous iterate.
C      45..47  
C          SIGH   INTERN Storage to hold raw error estimates of latest
C                        three iterates for mean error calculation.
C      48..50            Reserved
C      
C     Error codes (Stored to IERR):
C     =============================
C
C      -1     Iteration number limit (as set by LITMAX) exceeded 
C      -2     Iteration diverges: The error estimate exceeds an
C             internal maximum allowed value (DIFLIM).
C      -3     Iteration diverges: Bad tau values and norms of 
C             corrections oscillate in a two-iterates-cycle.
C      -4     Iteration cannot continue: bad tau value immediate 
C             at start or restart.
C      -10    Insufficient workspace  
C      -20    Bad input to dimensional parameter N
C      -21    Nonpositive value for RTOL supplied
C      -22    Negative scaling value via vector XSCAL supplied
C      -25    Block diagonal scaling selected, and N is not a multiple
C             of the specified diagonal blocks size NBSIZE.
C      -30    One or more fields specified in IOPT are invalid
C             (for more information, see error-printout)
C
C*    Machine dependent constants used:
C     =================================
C
C     GREAT = squareroot of maxreal divided by 10
C     SMALL = squareroot of "smallest positive machine number
C             divided by relative machine precision"
C     DOUBLE PRECISION GREAT,SMALL  (used in GBPCHK)
C     ------------------------------------------------------------
C*    End Prologue
      EXTERNAL MULJAC, PRECON
      INTEGER MPRMON, LUMON, NBSIZE
C
      LUMON = IOPT(14)
      IF (LUMON.LT.1 .OR. LUMON.GT.99) LUMON=6
      MPRMON = IOPT(13)
C
      NZMAX = INDAWK(1)
      IF (NZMAX.NE.0) THEN
        IF (NZMAX.LT.NELT) THEN
          IERR = -26
          IF (MPRMON.GT.0)  WRITE(LUMON, 10001) NZMAX, NELT
10001     FORMAT(1X,'Invalid NZMAX (=INDAWK(1)) specified: ',I7,
     $           ' - must be greater than NELT = ',I8,/
     $           ' or zero if equal to NELT')
          RETURN
        ELSE IF (NZMAX.GT.NELT) THEN 
          DO 10 I=1,NELT
            INDAWK(11+NELT+I) = INDAWK(11+NZMAX+I) 
10        CONTINUE
        ENDIF
      ENDIF
      INDAWK(1) = LA
      INDAWK(2) = LINDAW
      INDAWK(3) = 0
      INDAWK(4) = 1
      INDAWK(5) = NELT
      INDAWK(7) = LUMON
      INDAWK(8) = MPRMON
      INDAWK(9) = IOPT(35)
      IF ( IOPT(35) .EQ. 4 ) THEN
        NBSIZE = IOPT(36)
        IF (NBSIZE.EQ.0) NBSIZE = 2
        IF ( MOD(N,NBSIZE) .NE. 0 ) THEN
          IERR = -25
          IF (MPRMON.GT.0)  WRITE(LUMON, 10002) N, NBSIZE
10002     FORMAT(1X,'Incompatible size N of the linear system and ',
     $            /,'block size NBSIZE for block diagonal scaling',/,
     $              'N must be a multiple of NBSIZE  -  ',
     $              'N = ',I7,'   NBSIZE = ',I4)
          RETURN 
        ENDIF
        INDAWK(10) = NBSIZE
      ENDIF
      INDAWK(11) = NELT
C     Setup preconditioner by EASYPACK subroutine EJAC
      CALL EJAC(DUMMY,N,X,XSCAL,RDUMMY,A,INDAWK,IDUMMY,IERR)
      IF (IERR.NE.0) RETURN
C     Call Good Broyden
      IOPT(35) = 0
      IOPT(36) = 0
      CALL GBIT1(N,MULJAC,PRECON,RHS,X,XSCAL,RTOL,IOPT,IERR,
     $                 LRWK,RWK,LINDAW,INDAWK,LA,A)
      RETURN
      END
C
C
      SUBROUTINE JAC1(FCN,N,U,UWGT,F,NELT,IDUMMY,A,IA,JA,NFILL,
     $                RWKU,IWKU,NJAC,IFAIL)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      EXTERNAL FCN
      INTEGER N
      DOUBLE PRECISION U(N),UWGT(N),F(N)
      INTEGER NELT,IDUMMY
      DOUBLE PRECISION A(NELT)
      INTEGER IA(NELT),JA(NELT)
      INTEGER NFILL
      DOUBLE PRECISION RWKU(*)
      INTEGER IWKU(*)
      INTEGER NJAC,IFAIL 
C
C     Pass number of nonzero matrix elements to EASYPACK
      NFILL = IWKU(1)
      RETURN
      END
C
C
      SUBROUTINE FCN1()
C
C     Dummy external to satisfy reference from EASYPACK (not called)
C
      RETURN
      END
