      SUBROUTINE GBIT1(N,MULJAC,PRECON,RHS,X,XSCAL,RTOL,IOPT,IERR,
     $                 LRWK,RWK,LIWKU,IWKU,LRWKU,RWKU)
C*    Begin Prologue GBIT1
      INTEGER N
      EXTERNAL MULJAC,PRECON
      DOUBLE PRECISION RHS(N),X(N),XSCAL(N),RTOL
      INTEGER IOPT(50),IERR
      INTEGER LRWK
      DOUBLE PRECISION RWK(LRWK)
      INTEGER LIWKU
      INTEGER IWKU(LIWKU)
      INTEGER LRWKU
      DOUBLE PRECISION RWKU(LRWKU)
C     ------------------------------------------------------------
C
C*  Title
C
C     Good Broyden -
C     Iterative solution of a linear system (implicitly given matrix
C     by subroutine MULJAC - computes A*vector)
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*    External subroutines (to be supplied by the user)
C     =================================================
C
C     (Caution: Arguments declared as (input) must not
C               be altered by the user subroutines ! )
C
C       MULJAC( N, X, Y, RWKU, IWKU )
C                Ext    Jacobian * vector product subroutine
C         N        Int    The number of vector components (input)
C         X(N)     Dble   The vector to be multiplied by the Jacobian
C                         (input)
C         Y(N)     Dble   The array to get the result vector
C                         Jacobian * X (output)
C       * RWKU(*)  Dble   Real Workspace for the user routines
C       * IWKU(*)  Int    Integer Workspace for the user routines
C         
C         
C       PRECON( N, R, Z, RWKU, IWKU )
C                Ext    Preconditioning solver subroutine for system
C                       M * Z = R
C         N        Int    The number of vector components (input)
C         R(N)     Dble   The right hand side of the system (input)
C         Z(N)     Dble   The array to get the solution vector (output)
C       * RWKU(*)  Dble   Real Workspace for the user routines
C       * IWKU(*)  Int    Integer Workspace for the user routines
C
C*    Input parameters (* marks inout parameters)
C     ===========================================
C
C     N          Int     Dimension of linear system
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     LIWKU          Int    Declared dimension of user function 
C                           integer workspace.
C     IWKU(LIWKU)    Int    User integer workspace (passed to
C                           user subroutines MULJAC and PRECON).
C     LRWKU          Int    Declared dimension of user function 
C                           real workspace.
C     RWKU(LRWKU)    Dble   User real workspace (passed to
C                           user subroutines MULJAC and PRECON).
C
C     Note 1.
C        The iteration stops, if the criterion 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                          by the iterative linear solver Good Broyden.
C                          Values <=0 will be special handled as listed
C                          below:
C                          a. An input <=-2 means KMAX=0 will be used. 
C                          b. A "-1" input means, that there no limit 
C                             applies, e.g. the Good Broyden will be
C                             served with the maximum possible value 
C                             allowed by the total workspace amount (but
C                             a value < 2 will not be accepted!)
C                             On output, KMAX will be set to the compu-
C                             ted value, which has (or would have) been
C                             accepted
C                          c. A zero input means KMAX will be set to the
C                             default value 9.
C      32  LITMAX 100      Maximum number of iterations allowed
C      33..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*    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      -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
      INTEGER L
      PARAMETER ( L=3 )
C     (Internal note: For L>6, RWK-Base needs to be increased!)
      INTEGER KMAXDF,KMAXMI
      PARAMETER ( KMAXDF=9, KMAXMI=2 )
      DOUBLE PRECISION ZERO,RHOS,TAUMIN,TAUMAX,TAUEQU
      PARAMETER( ZERO=0.0D0 )
      PARAMETER( RHOS=4.0D0 )
      PARAMETER( TAUMIN=1.0D-8, TAUMAX=1.0D2 , TAUEQU=1.0D-2 )
      LOGICAL QSUCC,QINIMO
      INTEGER LUMON, MPRMON, KMAX, IDMUL, IDPRE,
     $        L1, L2, L4, L5, L6, L7, L8
      DOUBLE PRECISION DELNRM
      EXTERNAL GBIT1I
      CHARACTER CHGDAT*20, PRODCT*8
C
C     Version: 1.0               Latest change:
C     -----------------------------------------
C
      DATA      CHGDAT      /'March 27, 1991      '/
      DATA      PRODCT      /'GBIT1   '/
      IERR = 0
      QSUCC = IOPT(1).EQ.1
C        Print iteration monitor?
      MPRMON= IOPT(13)
      LUMON= IOPT(14)
      IF (LUMON .LE. 0 .OR. LUMON .GT. 99) THEN
        LUMON = 6
        IOPT(14)=LUMON
      ENDIF
      QINIMO = MPRMON.NE.0.AND..NOT.QSUCC
C     Print GBIT1 heading lines
      IF(QINIMO)THEN
10000   FORMAT('   G B I T 1 *****  V e r s i o n  ',
     $         '1 . 1 ***',//,1X,'Good Broyden ',
     $         'for the iterative solution of linear systems',//)
        WRITE(LUMON,10000)
      ENDIF
C     Check input parameters and options
      CALL GBPCHK(N,X,XSCAL,RTOL,IOPT,IERR,LRWK,RWK)
C     Exit, if any parameter error was detected till here
      IF (IERR.NE.0) RETURN 
C        Print time summary statistics?
      MPRTIM = IOPT(19)
      LUTIM = IOPT(20)
C
      IF (IOPT(31).EQ.0) IOPT(31)=KMAXDF
      KMAX = IOPT(31)
      IF (IOPT(31).LT.-1) KMAX = 0
      IF (IOPT(31).EQ.-1) THEN
        KMAX = IDINT( DBLE(FLOAT(LRWK-5*N-50)) / DBLE(FLOAT(N+2)))
        KMAX = MAX0(KMAX,0)
        IF ( KMAX.LT.KMAXMI ) THEN
          IF (MPRERR.GE.1) WRITE (LUERR,10005) KMAX
10005     FORMAT(1X,'Workspace optimal KMAX would be ',I7,
     $              ' - but is too small')
          KMAX = KMAXMI
        ENDIF 
        IOPT(31) = KMAX
        IF (MPRMON.GE.1) WRITE (LUMON,10006) KMAX 
10006   FORMAT(1X,'selected KMAX is ',I7)
      ENDIF
      L1 = 51
      L2 = L1 + N
      L3 = L2 + N*KMAX
      L4 = L3 + N
      L5 = L4 + N
      L6 = L5 + N
      L61 = L6 + N
      L7 = L61 + KMAX
      L8 = L7 + KMAX
      NRW = L8 - 1
      IOPT(44)=NRW
      IOPT(45)=L2-1
      IF (MPRMON.LT.0) THEN
        WRITE(LUMON,1000) KMAX 
1000    FORMAT(6X,'KMax = ',I5)     
      ENDIF
      IF (NRW.GT.LRWK) THEN
        IERR = -10
        IF (MPRMON.NE.0) WRITE(LUMON,10020) 'Real',NRW-LRWK
10020   FORMAT(' GBIT1 - ',A,' Workspace exhausted,',
     $         ' at least more required: ',I6)
        RETURN
      ENDIF
C
      IF (IOPT(21).EQ.0) THEN
        IOPT(22)=0
      ELSE
        IF (IOPT(22).EQ.0) IOPT(22)=6
      ENDIF
C
      LITMAX=IOPT(32)
      IF (LITMAX.EQ.0) LITMAX=100
      IOPT(32)=LITMAX
C     Use always standard termination criterion
      IOPT(33)=1
C
      IF (RWK(21).EQ.ZERO) RWK(21)=RHOS
      IF (RWK(22).EQ.ZERO) RWK(22)=TAUMIN
      IF (RWK(23).EQ.ZERO) RWK(23)=TAUMAX
      IF (RWK(24).EQ.ZERO) RWK(24)=TAUEQU
C
      IF (.NOT.QSUCC) THEN
        DO 20 J=1,N
          RWK(L1-1+J)=ZERO
20      CONTINUE
      ENDIF
C
C     Initialize and start time measurements monitor
C
      IF ( MPRTIM.NE.0 ) THEN
        CALL MONINI (' GBIT1',LUTIM)
        CALL MONDEF (0,'GBIT1')
        CALL MONDEF (1,'Muljac')
        CALL MONDEF (2,'Precon')
        CALL MONSTR (IERR)
      ENDIF
      IDMUL = 1
      IDPRE = 2
C
      CALL GBIT1I(N,LITMAX,RTOL,XSCAL,MULJAC,PRECON,
     $           X,RWK(L1),RHS,IOPT,QSUCC,IERR,
     $           KMAX,RWK(L2),RWK(L3),RWK(L4),RWK(L5),
     $           RWK(L6),RWK(L61),RWK(L7),IOPT(41),IOPT(42),
     $           L,RWK(45),RWK(41),RWK(42),RWK(43),RWK(44),
     $           LIWKU,IWKU,LRWKU,RWKU,DELNRM,RWK(32),RWK(31),
     $           RWK(21),RWK(22),RWK(23),RWK(24),IOPT(43),IDMUL,IDPRE)
C
      IF ( MPRTIM.NE.0 ) THEN
         CALL MONHLT
         CALL MONPRT
       ENDIF

C
C     End of subroutine GBIT1
      RETURN
      END
C
      SUBROUTINE GBPCHK(N,X,XSCAL,RTOL,IOPT,IERR,LRWK,RWK)
C*    Begin Prologue GBPCHK
      INTEGER N
      DOUBLE PRECISION X(N),XSCAL(N)
      DOUBLE PRECISION RTOL
      INTEGER IOPT(50)
      INTEGER IERR
      INTEGER LRWK
      DOUBLE PRECISION RWK(LRWK)
C     ------------------------------------------------------------
C
C*    Summary :
C
C     G B P C H K : Checking of input parameters and options
C                   for GBIT1.
C
C*    Parameters:
C     ===========
C
C     See parameter description in driver routine.
C
C*    Subroutines called: ZIBCONST
C
C*    Machine dependent constants used:
C     =================================
C
C     EPMACH = relative machine precision
C     GREAT = squareroot of maxreal divided by 10
C     SMALL = squareroot of "smallest positive machine number
C             divided by relative machine precision"
      DOUBLE PRECISION EPMACH,GREAT,SMALL
C
C     ------------------------------------------------------------
C*    End Prologue
C
      INTRINSIC DBLE
      DOUBLE PRECISION ONE
      PARAMETER (ONE=1.0D0)
      DOUBLE PRECISION TEN
      PARAMETER (TEN=1.0D1)
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
C
      PARAMETER (NUMOPT=50)
      INTEGER IOPTL(NUMOPT),IOPTU(NUMOPT)
      DOUBLE PRECISION TOLMIN,TOLMAX,DEFSCL
C
      DATA IOPTL /0,0,0,0,0,0,0,0,0,0,0,0,-9999999,1,0,0,0,0,0,1,
     $            0,1,0,0,0,0,0,0,0,0,-2,0,0,0,0,
     $            0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/
      DATA IOPTU /1,0,0,0,0,0,0,0,0,0,0,0,1,99,0,0,1,0,1,99,
     $            1,99,0,0,0,0,0,0,0,0,9999999,9999999,0,0,0,
     $            0,0,0,0,0,9999999,9999999,9999999,9999999,9999999,
     $            0,0,0,0,0/
C
      CALL ZIBCONST(EPMACH,SMALL)
      GREAT  = TEN/SMALL
      IERR = 0
C        Print error messages?
      MPRERR = IOPT(13)
      LUERR = IOPT(14)
C     Default output units
      IF (LUERR .LE. 0 .OR. LUERR .GT. 99) THEN
        LUERR = 6
        IOPT(14)=LUERR
      ENDIF
      IF (IOPT(20).EQ.0) IOPT(20)=6
      IF (IOPT(22).EQ.0) IOPT(22)=6
C
C     Checking dimensional parameter N
      IF ( N.LE.0 ) THEN
        IF (MPRERR.NE.0)  WRITE(LUERR,10011) N
10011   FORMAT(/,' Error: Bad input to dimensional parameter N supplied'
     $         ,/,8X,'choose N positive, your input is: N = ',I5)
        IERR = -20
      ENDIF
C
C     Checking and conditional adaption of the user-prescribed RTOL
      IF (RTOL.LE.ZERO) THEN
        IF (MPRERR.NE.0) 
     $      WRITE(LUERR,'(/,A)') ' Error: Nonpositive RTOL supplied'
        IERR = -21
      ELSE
        TOLMIN = EPMACH*TEN*DBLE(N)
        IF(RTOL.LT.TOLMIN) THEN
          RTOL = TOLMIN
          IF (MPRERR.LT.0) 
     $      WRITE(LUERR,10012) 'increased ','smallest',RTOL
        ENDIF
        TOLMAX = 1.0D0
        IF(RTOL.GT.TOLMAX) THEN
          RTOL = TOLMAX
          IF (MPRERR.LT.0) 
     $      WRITE(LUERR,10012) 'decreased ','largest',RTOL
        ENDIF
10012   FORMAT(/,' Warning: User prescribed RTOL ',A,'to ',
     $         'reasonable ',A,' value RTOL = ',D11.2)
      ENDIF
C     
C     Test user prescribed accuracy and scaling on proper values
      IF (N.LE.0) RETURN 
      DEFSCL = ONE
      DO 10 I=1,N
        IF (XSCAL(I).LT.ZERO) THEN
          IF (MPRERR.NE.0) THEN 
            WRITE(LUERR,10013) I
10013       FORMAT(/,' Error: Negative value in XSCAL(',I5,') supplied')
          ENDIF
          IERR = -22
        ENDIF
        IF (XSCAL(I).EQ.ZERO) XSCAL(I) = DEFSCL
        IF ( XSCAL(I).GT.ZERO .AND. XSCAL(I).LT.SMALL ) THEN
          IF (MPRERR.LT.0) THEN
            WRITE(LUERR,10014) I,XSCAL(I),SMALL
10014       FORMAT(/,' Warning: XSCAL(',I5,') = ',D9.2,' too small, ',
     $             'increased to',D9.2)
          ENDIF
          XSCAL(I) = SMALL
        ENDIF
        IF (XSCAL(I).GT.GREAT) THEN
          IF (MPRERR.LT.0) THEN
            WRITE(LUERR,10015) I,XSCAL(I),GREAT
10015       FORMAT(/,' Warning: XSCAL(',I5,') = ',D9.2,' too big, ',
     $             'decreased to',D9.2)
          ENDIF
          XSCAL(I) = GREAT
        ENDIF
10    CONTINUE
C     Checks options
      DO 20 I=1,NUMOPT
        IF (IOPT(I).LT.IOPTL(I) .OR. IOPT(I).GT.IOPTU(I)) THEN
          IERR=-30
          IF (MPRERR.NE.0) THEN
            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)
          ENDIF
        ENDIF
20    CONTINUE
C     End of subroutine GBPCHK
      RETURN
      END
C
      SUBROUTINE GBIT1I(N,LITMAX,RTOL,XSCAL,MULJAC,PRECON,
     $                 X,DEL,RHS,IOPT,QSUCC,IERR,KMAX,DELTA,DELTAK,
     $                 Q,Z,SITER,T,SIGMA,
     $                 ITER,K,L,SIGH,SIGMAK,TAUN1,TAUN2,SOLNRA,
     $                 LIWKU,IWKU,LRWKU,RWKU,
     $                 DELNRM,SOLNRM,DIFNRM,RHOS,TAUMIN,TAUMAX,TAUEQU,
     $                 NMULJ, IDMUL, IDPRE)
C*    Begin Prologue GBIT1I
      INTEGER N,LITMAX
      DOUBLE PRECISION RTOL,XSCAL(N)
      EXTERNAL MULJAC,PRECON
      INTEGER IOPT(50)
      LOGICAL QSUCC
      INTEGER IERR,KMAX
      DOUBLE PRECISION X(N),DEL(N),RHS(N),DELTA(N,KMAX),DELTAK(N),
     $                 Q(N),Z(N),SITER(N),T(KMAX),
     $                 SIGMA(KMAX)
      INTEGER ITER,K
      DOUBLE PRECISION SIGH(L),SIGMAK,TAUN1,TAUN2,SOLNRA
      INTEGER LIWKU
      INTEGER IWKU(LIWKU)
      INTEGER LRWKU
      DOUBLE PRECISION RWKU(LRWKU)
      DOUBLE PRECISION DELNRM,SOLNRM,DIFNRM,RHOS,TAUMIN,TAUMAX,TAUEQU
      INTEGER NMULJ,IDMUL,IDPRE
C     ------------------------------------------------------------
C
C*    Internal core routine
C
C     Iterative solution of a linear system (implicit given matrix
C     by subroutine MULJAC - computes A*vector)
C
C     ------------------------------------------------------------
C*    End Prologue
      DOUBLE PRECISION ZERO,ONE,GREAT
      INTEGER L
      PARAMETER( ZERO=0.0D0, ONE = 1.0D0, GREAT=1.0D+35 )
      EXTERNAL SPRODS
      INTRINSIC DABS,DSQRT
      INTEGER I,J,K1,IFINI,MPRMON,ISYM,LUMON,IPLOT,ITERM
      DOUBLE PRECISION DN,EPMACH,SMALL,FAKTOR,FAKT1,DIFLIM,
     $                 SPRODS,SIGMAP,GAMMA,TAU,TK,
     $                 SIGQN,DIFQN,
     $                 TKRMI,TKRMA
      DOUBLE PRECISION SIGNEW,SIGQ1,SIGMIN
      LOGICAL QSTOP,QSTOPK
C*    Begin
C
C  Initiation
C------------
C
      CALL ZIBCONST(EPMACH,SMALL)
      IERR=10000
      ISYM=0
      ITERM=0
C
      IFINI =IOPT(33)
      MPRMON=IOPT(13)
      LUMON =IOPT(14)
      IPROPT=IOPT(17)
      IPLOT =IOPT(22)
      MPRTIM=IOPT(19)
      TKRMI=ONE
      TKRMA=ONE
      SIGQN=GREAT
      DIFQN=GREAT
      DIFNRM=GREAT
      SOLNRM=GREAT
      TAU=GREAT
      DIFLIM=1.0D10
      QSTOP = .FALSE.
      QSTOPK = .FALSE.
      DN = DBLE(FLOAT(N))
      SIGMIN=DN*EPMACH**2
C
C  initial preparations
C----------------------
C
      DO 2 I=1,N
        SITER(I) = X(I) / XSCAL(I)
        X(I) = (X(I)-DEL(I)) / XSCAL(I)
        DEL(I) =DEL(I) / XSCAL(I)
2     CONTINUE
C
C     continuation entry
C
      IF (QSUCC) THEN
        SOLNRM=SOLNRA
        GOTO 20
      ENDIF
C
C   initiation for new iteration run
C   --------------------------------
C
      TAUN1=GREAT
      TAUN2=GREAT
      NMULJ=0
C
C  initial print
C
      IF (MPRMON.LT.0) THEN
        SOLNRM = DSQRT( SPRODS(N, X, X)/DN )
        IF (IPROPT.EQ.0) WRITE(LUMON,10020)
10020   FORMAT(2X,' It',7X,'Cor',7X,'Sol',1X,'EstAbsErr',6X,'SolQ',
     $         7X,'tau')
        IF (IPROPT.EQ.1) WRITE(LUMON,10030)
10030   FORMAT(2X,' It',7X,'Cor',7X,'Del',7X,'Sol',1X,'EstAbsErr',
     $           6X,'DelQ',6X,'SolQ',7X,'tau')
      ENDIF
      ITER = 0
      IOPT(49) = IOPT(49) + 1
C
C  start / entry for restart of iteration
C----------------------------------------
C
3     CONTINUE
C
C     k := 0 
      K = 0
C
C------------------------------------------------------------
C     --- r0 := b-A*x0   ===  z := rhs-A*siter ---
C     --- delta0 := H0*r0  ===  Solve A(precon)*deltak = z ---
C     --- sigma0 := delta0(T)*delta0 === sigmak := deltak(T)*deltak ---
C------------------------------------------------------------
C
      DO 8 I=1,N
        Z(I) = SITER(I)*XSCAL(I)
8     CONTINUE
      IF (MPRTIM.NE.0) CALL MONON (IDMUL)
      CALL MULJAC(N, Z, Q, RWKU, IWKU)
      IF (MPRTIM.NE.0) CALL MONOFF (IDMUL)
C
      DO 10 I=1,N
        Q(I) = RHS(I) - Q(I)
10    CONTINUE
C
      IF (MPRTIM.NE.0) CALL MONON (IDPRE)
      CALL PRECON(N, Q, DELTAK, RWKU, IWKU)
      IF (MPRTIM.NE.0) CALL MONOFF (IDPRE)
      NMULJ=NMULJ+1
      DO 15 I=1,N
        DELTAK(I) = DELTAK(I)/XSCAL(I)
15    CONTINUE
C 
      SIGMAK = SPRODS(N, DELTAK, DELTAK)
      IF (SIGMAK.LE.SIGMIN) GOTO 1000
C
      CALL ITZMID(L,ITER,SIGH,SIGMAK,SIGQN,.FALSE.)
C
C  Main iteration loop         
C---------------------
C
20    IF ( .NOT. QSTOPK .OR. .NOT. QSTOP ) THEN
C
        IF (ITER.GE.LITMAX) THEN
          IERR = -1
          GOTO 998 
        ENDIF
        K1=K+1
C
C------------------------------------------------------------
C    --- qk := A*deltak  ===  q := A*deltak ---
C    --- z0quer := H0*qk  ===  Solve A(precon)*z = q ---
C------------------------------------------------------------
C
        DO 30 I=1,N
          Z(I) = DELTAK(I)*XSCAL(I)
30      CONTINUE
        IF (MPRTIM.NE.0) CALL MONON (IDMUL)
        CALL MULJAC(N, Z, Q, RWKU, IWKU)
        IF (MPRTIM.NE.0) CALL MONOFF (IDMUL)
C
        IF (MPRTIM.NE.0) CALL MONON (IDPRE)
        CALL PRECON(N, Q, Z, RWKU, IWKU)
        IF (MPRTIM.NE.0) CALL MONOFF (IDPRE)
        NMULJ=NMULJ+1
        DO 40 I=1,N
          Z(I) = Z(I)/XSCAL(I)
40      CONTINUE
C
C  update loop
C------------- 
C
        DO 100 I=1,K-1
          FAKTOR = SPRODS(N, DELTA(1,I), Z) / SIGMA(I)
          FAKT1 = ONE-T(I)
          DO 110 J=1,N
            Z(J) = Z(J) + FAKTOR * (DELTA(J,I+1) - FAKT1*DELTA(J,I))
110       CONTINUE
100     CONTINUE
        IF (K.NE.0) THEN
          FAKTOR = SPRODS(N, DELTA(1,K), Z) / SIGMA(K)
          FAKT1 = ONE-T(K)
          DO 120 J=1,N
            Z(J) = Z(J) + FAKTOR * (DELTAK(J) - FAKT1*DELTA(J,K))
120       CONTINUE
        ENDIF
C
C----------------------------------------------------------
C    --- zk := zquerk ===  z now corresponds to zk ---
C    --- gammak := deltak(t)*zk ---
C    --- tauk := sigmak/gammak ---
C----------------------------------------------------------
C
        GAMMA = SPRODS(N, DELTAK, Z)
        IF (DABS(GAMMA).LE.SIGMIN) GOTO 1001
        IF (GAMMA.NE.ZERO) TAU = SIGMAK / GAMMA
        TK = TAU
C
C  check for restart condition
C-----------------------------
C
        IF (GAMMA.EQ.ZERO .OR. TAU.LT.TAUMIN .OR. TAU.GT.TAUMAX) THEN
          IF (K.EQ.1 .AND. GAMMA.NE.ZERO ) THEN
            IF ( DABS(TAUN2-TAU) .LT. TAUEQU*DABS(TAU) ) THEN
              IERR=-3
              GOTO 998
            ENDIF
            TAUN2 = TAUN1
            TAUN1 = TAU
          ENDIF
          IF(MPRMON.LT.0.AND.GAMMA.NE.ZERO)
     1       WRITE(LUMON,10003) TAU,ITER,K
10003     FORMAT(' >>> Restart required due to tauk = ',D10.3,
     $           '( iter=',I5,',k=',I3,' )')
          IF(MPRMON.LT.0.AND.GAMMA.EQ.ZERO) 
     $      WRITE(LUMON,*)' >>> Restart required due to gammak = 0.0d0'
          IF (K.EQ.0) THEN
            IF (GAMMA.EQ.ZERO) THEN
              IF (MPRMON.LT.0)  WRITE(LUMON,10013) 
10013         FORMAT(' >>> Termination - restart not possible')
              IERR = -4
              GOTO 998
            ENDIF
            IF (MPRMON.LT.0) WRITE(LUMON,10023) TAU,ITER,K
10023       FORMAT(' >>> Restart condition ignored - tau = ',D10.3,
     $             '( iter=',I5,',k=',I3,' )')
            IF (TAU.LT.TAUMIN) TK = TKRMI
            IF (TAU.GT.TAUMAX) TK = TKRMA
            IF (MPRMON.LT.0) WRITE(LUMON,10033) TK
10033       FORMAT (' >>> tk reset to ',D10.3)
          ELSE
            GOTO 9990
          ENDIF
        ENDIF
C
C----------------------------------------------------------
C    --- x(k+1) := xk + tk*deltak ---
C    --- delta(k+1) := deltak - tauk*zk (if tk=tauk) ---
C                        - or -
C    --- delta(k+1) := (1-tk+tauk)*deltak - tauk*zk (if tk<>tauk)
C    --- sigma(k+1) := delta(k+1)(t)*delta(k+1) ---
C----------------------------------------------------------
C
C  update iterate
C---------------- 
C
        DO 130 J=1,N
          DEL(J) = DEL(J) + TK*DELTAK(J)
          SITER(J) = X(J) + DEL(J)
130     CONTINUE
C
C  compute norms for converge check 
C
        SOLNRA = SOLNRM
        SOLNRM = DSQRT( SPRODS(N, SITER, SITER)/DN )
        DELNRM = DSQRT (SPRODS(N, DEL, DEL)/DN )
        DIFNRM = DABS(TK)*DSQRT(SIGMAK/DN)
C
C  save information to perform next update loop
C
        IF (K1.LE.KMAX) THEN
          DO 6 I=1,N
            DELTA(I,K1)=DELTAK(I)
6         CONTINUE
          SIGMA(K1) = SIGMAK
          T(K1) = TK
        ENDIF
C
C  new delta
C
         IF (TK.EQ.TAU) THEN
           DO 140 J=1,N
             DELTAK(J) = DELTAK(J) -  TAU * Z(J)
140        CONTINUE
         ELSE
           FAKTOR = ONE - TK + TAU
           DO 145 J=1,N
             DELTAK(J) = FAKTOR*DELTAK(J) -  TAU * Z(J)
145        CONTINUE
         ENDIF
C
C  new sigma
C
        SIGMAP = SIGMAK
        SIGMAK = SPRODS(N, DELTAK, DELTAK)
        SIGNEW = DSQRT(SIGMAK/DN)
C
Ctaum;  CALL ITZMID(L,ITER,SIGH,TK**2*SIGMAK,SIGQN,.TRUE.)
        CALL ITZMID(L,ITER,SIGH,SIGMAK,SIGQN,.TRUE.)
C
        SIGQ1=DSQRT(SIGQN/DN)
Ctaus;  SIGQ1=TK*DSQRT(SIGMAK/DN)
CSim;   SIGQ1=DSQRT(SIGMAK/DN)
C
C
C       write graphics data
C
        IF (IPLOT.GT.0) 
     $   WRITE(IPLOT,9510) ITER,K,DIFNRM,SIGNEW,SIGQ1,SIGQ1,SIGQ1,SOLNRM
9510    FORMAT(I4,I3,6(1PD12.3))
C
        DIFQN=RHOS*SIGQ1
C
C       --- Print monitor ---
C
        IF (MPRMON.LT.0) THEN
10004     FORMAT(1X,I4,7(1X,D9.2))
          IF (MOD(ITER,-MPRMON).EQ.0) THEN 
            IF (IPROPT.EQ.0)
     $          WRITE(LUMON,10004) ITER, DIFNRM, SOLNRA, DIFQN, SOLNRM,
     $                              TAU
            IF (IPROPT.EQ.1)
     $          WRITE(LUMON,10004) ITER, DIFNRM, DELNRM, SOLNRA,
     $                              DIFQN, DELNRM, SOLNRM, TAU 
          ENDIF 
        ENDIF
C
C  check for termination
C-----------------------
C  
        QSTOPK=QSTOP
C??!    QSTOPK=.TRUE.
        ITERM = 0
        IF (SIGMAK.LE.SIGMIN) QSTOPK=.TRUE. 
        IF (ITER.GE.L-1) THEN
CSim;   IF (ITER.GE.0) THEN
          IF (DIFQN.GT.GREAT*SOLNRM) IERR=-2
          IF (DIFQN.GT.DIFLIM) IERR=-2
          IF (IERR.EQ.-2) GOTO 998
          QSTOP = DIFQN .LE. RTOL*SOLNRM
          IF (QSTOP) ITERM = 1
          IF (IFINI.EQ.2)  QSTOP = DIFQN .LE. 0.25D0 * DELNRM .OR. QSTOP
          IF (ITERM.EQ.0 .AND. QSTOP) ITERM = 2
        ELSE
          QSTOP=.FALSE.
        ENDIF
        IF (SIGMAK.LE.SIGMIN) QSTOP=.TRUE. 
C 
C  decision for next step (proceed or restart due to kmax)
        K = K+1
        ITER = ITER+1
        IF (K.LE.KMAX) GOTO 20
C
C --- End of Iteration Loop ---
C
        IF (MPRMON.LT.0.AND.KMAX.GT.0) 
     $    WRITE(LUMON,10050) ITER
10050 FORMAT(' >>> Restart due to k > kmax ( iter=', I5, ' )')
C
C --- Restart ---
C
9990    CONTINUE
        IOPT(48) = IOPT(48) + (K*(K-5))/2 + 1 
        IOPT(49) = IOPT(49) + 1
        IF ( SIGMAK.GT.SIGMIN ) GOTO 3
      ENDIF
C
C  solution exit
C
1000  CONTINUE
      IF ( SIGMAK.GT.SIGMIN ) IOPT(48) = IOPT(48) + (K*(K-5))/2 + 1
      IERR = 0
      IF (MPRMON.LT.0)  WRITE(LUMON,*) ' > Solution found'
      GOTO 999
1001  CONTINUE
      IERR = 0
      IF (MPRMON.LT.0)  THEN
         WRITE(LUMON,*) ' > Solution found ???'
         WRITE(LUMON,*) ' sigmak/gamma:',SIGMAK,GAMMA
      ENDIF
      SIGMAK=GAMMA
      GOTO 999
C
C  emergency termination
C
998   CONTINUE
      IF (MPRMON.NE.0) THEN
        IF (IERR.EQ.-1) WRITE(LUMON,99801)  
99801   FORMAT(' >> Iter. number limit reached')
        IF (IERR.EQ.-2)  WRITE(LUMON,99802) 
99802   FORMAT(' >> Termination, since iteration diverges')
        IF (IERR.EQ.-3)  WRITE(LUMON,99803) TAUN1,TAU
99803   FORMAT(' >> Term., since bad tau oscillates, ',
     $         ' taui = ',D9.2,' , ',D9.2)
        IF (IERR.EQ.-4)  WRITE(LUMON,99804) 
99804   FORMAT(' >> Termination, since negative tau at restart')
      ENDIF
C
C  final update of return arguments
C
999   CONTINUE
      DO 9991 J=1,N
        X(J)=SITER(J)*XSCAL(J)
        DEL(J)=DEL(J)*XSCAL(J)
9991  CONTINUE
C
      IF (SIGMAK.LE.SIGMIN) THEN
        DIFNRM=SIGMAK
        RTOL=SIGMAK
        IOPT(50) = 0
      ELSE
        RTOL = DIFQN/SOLNRM
        DIFNRM = DABS(TAU)*DSQRT(SIGMAP/DN)
        IOPT(50) = ITERM
      ENDIF
      IF (ITER.NE.0) SOLNRA = SOLNRM
C
      IF (MPRMON.NE.0 .AND. IOPT(1).EQ.0)
     1   WRITE(LUMON,1090) ITER,RTOL,IOPT(50)
1090  FORMAT(10X,'LinSol:   ','Iter:',I6,'  EstPrec:',D11.2,
     1       '  TermCrit:',I2)
      IF (MPRMON.NE.0 .AND. IOPT(1).EQ.1)
     1   WRITE(LUMON,1091) ITER,RTOL,IOPT(50)
1091  FORMAT(10X,' Cont.:   ','Iter:',I6,'  EstPrec:',D11.2,
     1       '  TermCrit:',I2)
C
      RETURN
C
C  End of subroutine GBIT1I
C
      END
C
      SUBROUTINE ITZMID(L,ITER,DIFH,DIFNRM,DIFQN,QNEXT)
      INTEGER L,ITER
      DOUBLE PRECISION DIFH(L),DIFNRM,DIFQN
      LOGICAL QNEXT
      IF (ITER.LT.L) THEN
        DIFH(ITER+1) = DIFNRM
      ELSE
        IF (QNEXT) THEN
          DO 119 J=1,L-1
            DIFH(J)=DIFH(J+1)
119       CONTINUE
        ENDIF
        DIFH(L)=DIFNRM
      ENDIF
      IF (ITER.GE.L-1) THEN
        DIFQN = DIFH(1)+DIFH(L)
        DO 1191 J=2,L-1
          DIFQN = DIFQN+2.0D0*DIFH(J)
1191    CONTINUE
        DIFQN = DIFQN/DBLE(FLOAT(L+1))
      ENDIF
      RETURN
      END
C
      DOUBLE PRECISION FUNCTION SPRODS( N, X, Y )
      INTEGER N
      DOUBLE PRECISION X(N),Y(N)
      INTEGER I
      DOUBLE PRECISION S
C*    Begin
      S=0.0D0
      DO 10 I=1,N
        S = S + X(I)*Y(I)
10    CONTINUE
      SPRODS=S
      RETURN
      END
