PROGRAM TNLE1S
      IMPLICIT DOUBLEPRECISION(S)
C
C     ____________________________________________________________
C
C     Testexample for NLEQ1S: Chemical equilibrium problem.
C
C*  Written by        L. Weimann 
C*  Purpose           Testexample for code NLEQ1S
C*  Version           2.4
C*  Revision          March 2009
C*  Latest Change     March 2009
C*  Library           CodeLib
C*  Code              Fortran 77, Double Precision
C*  Environment       Standard Fortran 77 environment on PC's,
C                     workstations and hosts.
C
C     ____________________________________________________________
C
      INTEGER N, NFMAX, NBROY
      PARAMETER (N=6, NFMAX=18, NBROY=0)
      INTEGER LRWK
      PARAMETER (LRWK=6*NFMAX+(12+NBROY)*N+112)
      INTEGER LIWK
      PARAMETER (LIWK=12*NFMAX+11*N+193)
      INTEGER I
      DOUBLE PRECISION EPS
      INTEGER IOPT(50)
      INTEGER IERR,IFAIL
      DOUBLE PRECISION X(N),XSCAL(N),RW(LRWK)
      INTEGER IW(LIWK)
      REAL STIME,ETIME,CPTIME
      EXTERNAL F
      EXTERNAL DF
C:    Begin
      OPEN(2,FILE='nleq1s_analytic.dat')
      OPEN(9,FILE='nleq1s_analytic.out')
      WRITE(6,*) ' monitor: nleq1s_analytic.out,',
     $           ' data: nleq1s_analytic.dat'
1       FORMAT(' Nonlinear problems with sparse ',
     $  'Jacobian:',/,' Chemical equilibrium problem')
        WRITE(9,1)
        EPS = 1.0D-5
        DO 710 I=1,50
          IOPT(I)=0
710      CONTINUE
        DO 711 I=1,LIWK
          IW(I)=0
711     CONTINUE
        DO 713 I=1,LRWK
          RW(I)=0.0D0
713     CONTINUE
C       Execution mode: 0=Standard Mode, 1=Stepwise mode
        IOPT(2)=1
C       0 = Sparse pattern of Jacobian may vary 
C       1 = Fixed sparse pattern for all Jacobians to be computed
        IOPT(37)=0
C       Problem classification:
C       1 = linear , 2 = mildly nonlinear  3 = highly nonlinear
C       4 = extremely nonlinear
C       IOPT(31)=3
C       Automatic row scaling of linear system :
C       0 = allowed , 1 = inhibited
        IOPT(35)=0
C       Set MPRERR, MPRMON, MPRSOL, MPRTIM
        IOPT(11)=3
        IOPT(13)=3
        IOPT(15)=2
        IOPT(19)=1
C       Set print units LUERR, LUMON, LUSOL, LUTIM
        IOPT(12)=9
        IOPT(14)=9
        IOPT(16)=2
        IOPT(20)=9
C       Solution output format:
C       0=standard format, 1= GRAZIL readable output
        IOPT(46)=0
C       Override maximum allowed number of iterations:
C       IW(31)=200
C       Override starting damping factor:
C       RW(21)=1.0D0
C       Override minimal allowed damping factor:
C       RW(22)=1.0D-2
C       Starting vector for the iteration
        X(1)=3.0D-4
        X(2)=3.0D-4
        X(3)=27.5D0
        X(4)=3.0D-4
        X(5)=27.5D0
        X(6)=27.5D0
C        X(1)=1.0D0
C        X(2)=1.0D0
C        X(3)=1.0D0
C        X(4)=1.0D0
C        X(5)=1.0D0
C        X(6)=1.0D0
        DO 75 I=1,N
          XSCAL(I) = 0.0D0
75      CONTINUE
        IERR=-1
        I=0
        CALL ZIBSEC(STIME,IFAIL)
31      IF (IERR.EQ.-1) THEN
          CALL NLEQ1S(N,NFMAX,F,DF,X,XSCAL,EPS,IOPT,IERR,LIWK,IW,
     $    LRWK,RW)
C         Clear workspace declared not to be used
          NIFREE=IW(16)
          DO 311 K=NIFREE,LIWK
            IW(K)=0
311       CONTINUE
          NRFREE=IW(17)
          DO 313 K=NRFREE,LRWK
            RW(K)=0.0D0
313       CONTINUE
          I=I+1
32         FORMAT(' Returned from call ',I4,' of NLEQ1S')
          WRITE(9,32)I
C         IOPT(2)=0
          GOTO 31
        ENDIF
        CALL ZIBSEC(ETIME,IFAIL)
        CPTIME = ETIME-STIME
73      FORMAT(//,1X,'Time used =',F9.3,1X,'Sec',//,66('*'),
     *    /)
        WRITE(9,73)CPTIME
      END
      SUBROUTINE F(N,X,FX,IFAIL)
      IMPLICIT DOUBLEPRECISION(S)
      INTEGER N
      DOUBLE PRECISION X(N),FX(N)
      INTEGER IFAIL
C:    End Parameter
C:    Begin
      FX(1)=X(1)+X(2)+X(4)-1.0D-3
      FX(2)=X(5)+X(6)-55.0D0
      FX(3)=X(1)+X(2)+X(3)+2.0D0*X(5)+X(6)-110.001D0
      FX(4)=X(1)-1.0D-1*X(2)
      FX(5)=X(1)-1.0D4*X(3)*X(4)
      FX(6)=1.0D-14*X(5)-55.0D0*X(3)*X(6)
      RETURN
      END
      SUBROUTINE DF(N,X,DFX,IROW,ICOL,NFILL,IFAIL)
      IMPLICIT DOUBLEPRECISION(S)
      INTEGER N
      INTEGER NFILL
      DOUBLE PRECISION X(N)
      DOUBLE PRECISION DFX(NFILL)
      INTEGER  IROW(NFILL),ICOL(NFILL)
C:    End Parameter
      INTEGER NSPACE
C:    Begin
      NSPACE  = NFILL
      NFILL  = 18
      IF(NFILL.GT.NSPACE) RETURN
      DFX(1)=1.0D0
      IROW(1)=1
      ICOL(1)=1
      DFX(2)=1.0D0
      IROW(2)=1
      ICOL(2)=2
      DFX(3)=1.0D0
      IROW(3)=1
      ICOL(3)=4
      DFX(4)=1.0D0
      IROW(4)=2
      ICOL(4)=5
      DFX(5)=1.0D0
      IROW(5)=2
      ICOL(5)=6
      DFX(6)=1.0D0
      IROW(6)=3
      ICOL(6)=1
      DFX(7)=1.0D0
      IROW(7)=3
      ICOL(7)=2
      DFX(8)=1.0D0
      IROW(8)=3
      ICOL(8)=3
      DFX(9)=2.0D0
      IROW(9)=3
      ICOL(9)=5
      DFX(10)=1.0D0
      IROW(10)=3
      ICOL(10)=6
      DFX(11)=1.0D0
      IROW(11)=4
      ICOL(11)=1
      DFX(12)=-0.1D0
      IROW(12)=4
      ICOL(12)=2
      DFX(13)=1.0D0
      IROW(13)=5
      ICOL(13)=1
      DFX(14)=-1.0D4*X(4)
      IROW(14)=5
      ICOL(14)=3
      DFX(15)=-1.0D4*X(3)
      IROW(15)=5
      ICOL(15)=4
      DFX(16)=-55.0D0*X(6)
      IROW(16)=6
      ICOL(16)=3
      DFX(17)=1.0D-14
      IROW(17)=6
      ICOL(17)=5
      DFX(18)=-55.0D0*X(3)
      IROW(18)=6
      ICOL(18)=6
      RETURN
      END