C*********************************************
C                                            *
C   PROGRAM   : PCOMP                        * 
C   MODULE    : S (RUNTIME EXECUTIVE)        *
C   ABSTRACT  : FORTRAN PRECOMPILER          *
C   KEY WORD  : AUTOMATIC DIFFERENTIATION    *
C   SOURCE    : PCOMP 2.3 by M.LIEPELT       *
C               PCOMP 3.0 by M.DOBMANN       *
C   COPYRIGHT : C.TRASSL, K.SCHITTKOWSKI     *
C               MATHEMATISCHES INSTITUT,     *
C               UNIVERSITAET BAYREUTH,       *
C               D-95440 BAYREUTH, GERMANY    *
C   DATE      : JUNE 8, 1999                 *
C   VERSION   : 5.5                          *
C                                            *
C*********************************************
C
C
C
C     SUBROUTINE SYMPRP (SYMFIL,WA,LWA,IWA,LIWA,UWA,UIWA,IERR,MODE,
C    /                   NVAR,NFUNC)
C     INTEGER SYMFIL
C     INTEGER LWA,LIWA
C     DOUBLE PRECISION WA(LWA)
C     INTEGER IWA(LIWA)
C     INTEGER UWA,UIWA
C     INTEGER IERR,MODE,NVAR,NFUNC
C
C**********************************************************************
C
C   S Y M P R P   -   LOAD INTERMEDIATE CODE GENERATED BY SYMINP FROM
C                     SYMFIL INTO WORKING ARRAYS.
C
C   PARAMETERS:
C      SYMFIL    - INPUT DEVICE; THE INTERMEDIATE CODE GENERATED BY
C                  SYMINP WAS WRITTEN TO THIS FILE AND IS NOW LOADED.
C      WA(LWA)   - REAL WORKING ARRAY, REQUIRED BY SYMPRP. ON RETURN,
C                  WA() CONTAINS THE INTERMEDIATE CODE.
C      IWA(LIWA) - INTEGER WORKING ARRAY, CF. WA().
C      UWA,UIWA  - INDICATE THE ACTUAL SPACE OF WA() AND IWA() THAT
C                  HAS BEEN USED BY THE SUBROUTINE.
C      IERR      - THE PARAMETER SHOWS THE REASON FOR TERMINATING THE
C                  SUBROUTINE. ON RETURN IERR COULD CONTAIN THE FOLLOW-
C                  ING VALUES:
C                  IERR = 0 : SUCCESSFUL TERMINATION.
C                  IERR > 0 : AN ERROR HAS BEEN DETECTED. FOR FURTHER
C                             INFORMATION CF. SUBROUTINE SYMERR.
C      MODE      - THE PARAMETER IS USED FOR THE RESERVATION OF SPACE
C                  FOR THE HESSIAN MATRIX
C      NVAR      - ON RETURN, NVAR CONTAINS THE NUMBER OF VARIABLES ON
C                  FUNCTION INPUT FILE
C      NFUNC     - ON RETURN, NFUNC CONTAINS THE NUMBER OF FUNCTIONS ON
C                  INPUT FILE
C
C**********************************************************************
C
C     INTEGER I,PWA,PIWA,PX,PIX
C
C     INTEGER GSMDEP,GETIWA,HILF
C     PARAMETER (GSMDEP=10)
C
C     INTEGER HSMDEP
C     PARAMETER (HSMDEP=10)
C
C     INTEGER NOGRAD,GRAD,HESS
C     PARAMETER (NOGRAD=0,GRAD=1,HESS=2)
C
C     INTEGER IIS,VIS,IIC,VIC,IRC,VRC,IVA,VVA,IFN,XFN,VFN,VGR,VHE,VPF,IV
C     PARAMETER (IIS=1,VIS=2,IIC=3,VIC=4,IRC=5,VRC=6,IVA=7,VVA=8)
C     PARAMETER (IFN=9,XFN=10,VFN=11,VGR=12,VHE=13,VPF=14,IV=15)
C     INTEGER INFOLI(15)
C
C     DO 10 I=1,15
C       READ(SYMFIL,'(I6)',ERR=100) IWA(I)
C10   CONTINUE
C     PIWA=IWA(1)*5+IWA(2)+IWA(3)*4+IWA(4)+IWA(5)*4+IWA(7)*3+
C    1     IWA(9)*7+IWA(14)+15
C     PWA=IWA(6)
C     IF (MODE .EQ. NOGRAD) THEN
C       PX=IWA(8)+IWA(11)+GSMDEP+HSMDEP+2
C     ELSE IF (MODE .EQ. GRAD) THEN
C       PX=IWA(8)+IWA(11)+IWA(12)*IWA(8)+GSMDEP*IWA(8)+HSMDEP+1
C     ELSE IF (MODE .EQ. HESS) THEN
C       PX=IWA(8)+IWA(11)+IWA(12)*IWA(8)+IWA(13)*IWA(8)*IWA(8)+
C    1     GSMDEP*IWA(8)+HSMDEP*IWA(8)*IWA(8)
C     ELSE
C       IERR=61
C     ENDIF
C     IF (MODE .EQ. NOGRAD) IWA(12)=0
C     IF (MODE .NE. HESS) IWA(13)=0
C     PIX=IWA(15)
C     IF ((PWA+PX .GT. LWA) .OR. (PIWA+PIX .GT. LIWA)) THEN
C       IERR=32
C       RETURN
C     ENDIF
C     DO 20 I=16,PIWA
C       READ(SYMFIL,'(I6)',ERR=100) IWA(I)
C20   CONTINUE
C
C     NVAR=IWA(8)
C     NFUNC=IWA(10)
C     INFOLI(1)=15
C     INFOLI(2)=IWA(1)*5+INFOLI(1)
C     INFOLI(3)=IWA(2)+INFOLI(2)
C     INFOLI(4)=IWA(3)*4+INFOLI(3)
C     INFOLI(5)=IWA(4)+INFOLI(4)
C     INFOLI(7)=IWA(5)*4+INFOLI(5)
C     INFOLI(9)=IWA(7)*3+INFOLI(7)
C     INFOLI(14)=IWA(9)*7+INFOLI(9)
C     INFOLI(15)=IWA(14)+INFOLI(14)
C
C     INFOLI(6)=0
C     INFOLI(8)=IWA(6)+INFOLI(6)
C     INFOLI(11)=IWA(8)+INFOLI(8)
C     INFOLI(12)=IWA(11)+INFOLI(11)
C     INFOLI(13)=IWA(12)+INFOLI(12)
C
C     DO 50 I=1,IWA(9)
C       HILF=GETIWA(IFN,I,1,IWA,LIWA,INFOLI)
C       IF (HILF.EQ.1) THEN
C         NFUNC=NFUNC + GETIWA(IFN,I,3,IWA,LIWA,INFOLI)-1
C       ENDIF
C  50 CONTINUE
C
C     DO 30 I=1,PWA
C       READ(SYMFIL,'(D24.17)',ERR=100) WA(I)
C30   CONTINUE
C     IERR=0
C     UWA=PWA+PX
C     UIWA=PIWA+PIX
C     RETURN
C100  IERR=26
C     RETURN
C     END
C
C
C
      SUBROUTINE SYMFUN (X,N,F,M,ACTIVE,WA,LWA,IWA,LIWA,DFX,DFXLEN,IERR)
      INTEGER N,M
      DOUBLE PRECISION X(N),F(M)
      LOGICAL ACTIVE(M)
      INTEGER LWA,LIWA
      DOUBLE PRECISION WA(LWA)
      INTEGER IWA(LIWA)
      INTEGER DFXLEN
      INTEGER DFX(*)
      INTEGER IERR
C
C**********************************************************************
C
C   S Y M F U N   -   EVALUATE SYMBOLICALLY DEFINED FUNCTIONS.
C
C   PARAMETERS:
C      X(N)        - ON INPUT, THE ONE-DIMENSIONAL ARRAY X HAS TO
C                    CONTAIN THE ARGUMENT THE FUNCTIONS ARE TO BE
C                    COMPUTED AT.
C      F(M)        - ON RETURN, F CONTAINS THE VALUES OF THE ACTIVE
C                    FUNCTIONS AT ARGUMENT X.
C      ACTIVE(M)   - THE LOGICAL ARRAY SPECIFIES WHICH OF THE M 
C                    FUNCTIONS ARE TO BE COMPUTED (ACTIVE(K) = .TRUE.).
C      WA(LWA)     - REAL WORKING ARRAY, CONTAINS THE INTERMEDIATE CODE 
C                    GENERATED BY SYMINP.
C      IWA(LIWA)   - INTEGER WORKING ARRAY, CF. WA(LWA). 
C      DFX(DFXLEN) - THE ARRAY SPECIFIES WHICH FIRST AND SECOND
C                    DERIVATIVES ARE TO BE COMPUTED BY CONTAINING THE
C                    NUMBER OF THE VARIABLES
C      IERR        - THE PARAMETER SHOWS THE REASON FOR TERMINATING THE 
C                    SUBROUTINE. ON RETURN IERR COULD CONTAIN THE
C                    FOLLOWING VALUES:
C                    IERR = 0 : SUCCESSFUL TERMINATION.
C                    IERR > 0 : AN ERROR HAS BEEN DETECTED. FOR FURTHER
C                             INFORMATION CF. SUBROUTINE SYMERR.
C
C**********************************************************************
C
      INTEGER PIIS,PVIS,PIIC,PVIC,PIRC,PVRC,PIVA,PVVA,PIFN,PXFN,PVFN
C     INTEGER PVGR,PVHE,PVPF,PIV,LIIS,LVIS,LIIC,LVIC,LIRC,LVRC,LIVA
      INTEGER PVGR,PVHE,PVPF,LIIS,LVIS,LIIC,LVIC,LIRC,LVRC,LIVA
      INTEGER MPIIS,MPVIS,MPIIC,MPVIC,MPIRC,MPVRC,MPIVA,MPVVA,MPIFN
      INTEGER MPVFN,MPVPF,MPIV,MODE
      INTEGER LVVA,LIFN,LVFN,LVGR,LVHE,LVPF,LIV,LGST,LHST
C
      INTEGER IVAL
      DOUBLE PRECISION FVAL
      INTEGER I,J,K,L,PC
C
      INTEGER GSMDEP
      PARAMETER (GSMDEP=10)
C
      INTEGER NOGRAD
      PARAMETER (NOGRAD=0)
C
      IERR=0
C
      PIIS=IWA(1)
      PVIS=IWA(2)
      PIIC=IWA(3)
      PVIC=IWA(4)
      PIRC=IWA(5)
      PVRC=IWA(6)
      PIVA=IWA(7)
      PVVA=IWA(8)
      PIFN=IWA(9)
      PXFN=IWA(10)
      PVFN=IWA(11)
      PVGR=IWA(12)
      PVHE=1
      PVPF=IWA(14)
C     PIV=IWA(15)
C
      LIIS=16
      LVIS=LIIS+PIIS*5
      LIIC=LVIS+PVIS
      LVIC=LIIC+PIIC*4
      LIRC=LVIC+PVIC
      LVRC=1
      LIVA=LIRC+PIRC*4
      LVVA=LVRC+PVRC
      LIFN=LIVA+PIVA*3
      LVFN=LVVA+PVVA
      LVGR=LVFN+PVFN
      LVHE=LVGR+PVGR*PVVA
      LVPF=LIFN+PIFN*7
      LIV=LVPF+PVPF
      LGST=LVHE+PVHE*1
      LHST=LGST+GSMDEP*PVVA
C
      IF (N .NE. PVVA) THEN
        IERR=43
        RETURN
      ENDIF
      IF (M .NE. PVFN-(PIFN-PXFN)) THEN
        IERR=44
        RETURN
      ENDIF
      DO 10 I=1,N
        WA(LVVA+I-1)=X(I)
 10   CONTINUE
      K=0
      DO 20 I=1,PXFN
        IF (IWA(LIFN+(I-1)+(1-1)*PIFN) .EQ. 0) THEN
          K=K+1
          IF (ACTIVE(K)) THEN
            PC=IWA(LIFN+(I-1)+(6-1)*PIFN)
            MPIIS=MAX(1,IWA(1))
            MPVIS=MAX(1,IWA(2))
            MPIIC=MAX(1,IWA(3))
            MPVIC=MAX(1,IWA(4))
            MPIRC=MAX(1,IWA(5))
            MPVRC=MAX(1,IWA(6))
            MPIVA=MAX(1,IWA(7))
            MPVVA=MAX(1,IWA(8))
            MPIFN=MAX(1,IWA(9))
            MPVFN=MAX(1,IWA(11))
            MPVPF=MAX(1,IWA(14))
            MPIV =MAX(1,IWA(15))
            MODE=NOGRAD
            DO 30 L=LIV,LIV+MPIV-1
              IWA(L)=0
 30         CONTINUE
            CALL EVAL(PC,IVAL,FVAL,MODE,PVVA,1,
     1                MPIIS,MPVIS,MPIIC,MPVIC,MPIRC,MPVRC,MPIVA,MPVVA,
     2                MPIFN,MPVFN,PVHE,MPVPF,MPIV,IWA(LIIS),IWA(LVIS),
     3                IWA(LIIC),IWA(LVIC),IWA(LIRC),WA(LVRC),IWA(LIVA),
     4                WA(LVVA),IWA(LIFN),WA(LVFN),WA(LVGR),WA(LVHE),
     5                IWA(LVPF),IWA(LIV),WA(LGST),WA(LHST),DFX,DFXLEN,
     6                IERR)
            IF (IERR .NE. 0) RETURN
            F(K)=WA(LVFN+IWA(LIFN+(I-1)+(4-1)*PIFN)-1)
          ENDIF
        ELSE IF (IWA(LIFN+(I-1)+(1-1)*PIFN) .EQ. 1) THEN
          DO 15 J=1,IWA(LIFN+(I-1)+(3-1)*PIFN)
            K=K+1
            IF (ACTIVE(K)) THEN
              PC=IWA(LIFN+(I-1)+(6-1)*PIFN)
              MPIIS=MAX(1,IWA(1))
              MPVIS=MAX(1,IWA(2))
              MPIIC=MAX(1,IWA(3))
              MPVIC=MAX(1,IWA(4))
              MPIRC=MAX(1,IWA(5))
              MPVRC=MAX(1,IWA(6))
              MPIVA=MAX(1,IWA(7))
              MPVVA=MAX(1,IWA(8))
              MPIFN=MAX(1,IWA(9))
              MPVFN=MAX(1,IWA(11))
              MPVPF=MAX(1,IWA(14))
              MPIV =MAX(1,IWA(15))
              MODE=NOGRAD
              DO 40 L=LIV,LIV+MPIV-1
                IWA(L)=0
 40           CONTINUE
              IWA(LIV+IWA(LIFN+(I-1)+(2-1)*PIFN)-1)=J
              CALL EVAL(PC,IVAL,FVAL,MODE,PVVA,1,
     1                MPIIS,MPVIS,MPIIC,MPVIC,MPIRC,MPVRC,MPIVA,MPVVA,
     2                MPIFN,MPVFN,PVHE,MPVPF,MPIV,IWA(LIIS),IWA(LVIS),
     3                IWA(LIIC),IWA(LVIC),IWA(LIRC),WA(LVRC),IWA(LIVA),
     4                WA(LVVA),IWA(LIFN),WA(LVFN),WA(LVGR),WA(LVHE),
     5                IWA(LVPF),IWA(LIV),WA(LGST),WA(LHST),DFX,DFXLEN,
     6                IERR)
              IF (IERR .NE. 0) RETURN
              F(K)=WA(LVFN+IWA(LIFN+(I-1)+(4-1)*PIFN)-1+(J-1))
            ENDIF
 15       CONTINUE
        ENDIF
 20   CONTINUE
      RETURN
      END
C
C
C
      SUBROUTINE SYMGRA (X,N,F,M,DF,MMAX,ACTIVE,WA,LWA,IWA,LIWA,DFX,
     1                   DFXLEN,IERR)
      INTEGER N,M,MMAX
      DOUBLE PRECISION X(N),F(M),DF(MMAX,N)
      LOGICAL ACTIVE(M)
      INTEGER LWA,LIWA
      DOUBLE PRECISION WA(LWA)
      INTEGER IWA(LIWA)
      INTEGER DFXLEN
      INTEGER DFX(*)
      INTEGER IERR
C
C***********************************************************************
C
C   S Y M G R A   -   EVALUATE SYMBOLICALLY DEFINED FUNCTIONS AND
C                     CORRESPONDING GRADIENTS.
C
C   PARAMETERS:
C      X(N)       - ON INPUT, THE ONE-DIMENSIONAL ARRAY X HAS TO CONTAIN
C                   THE ARGUMENT THE FUNCTIONS ARE TO BE COMPUTED AT.
C      F(M)       - ON RETURN, F CONTAINS THE VALUES OF THE ACTIVE
C                   FUNCTIONS AT ARGUMENT X.
C      DF(MMAX,N) - ON RETURN, DF CONTAINS THE GRADIENTS OF THE ACTIVE
C                   FUNCTIONS AT ARGUMENT X. IN THE DRIVING PROGRAM, THE
C                   ROW DIMENSION OF DF HAS TO BE EQUAL TO MMAX.
C      ACTIVE(M)  - THE LOGICAL ARRAY SPECIFIES WHICH OF THE M FUNCTIONS
C                   ARE TO BE COMPUTED ( ACTIVE(K) = .TRUE. ).
C      WA(LWA)    - REAL WORKING ARRAY, CONTAINS THE INTERMEDIATE CODE
C                   GENERATED BY SYMINP.
C      IWA(LIWA)  - INTEGER WORKING ARRAY, CF. WA(LWA).
C      DFX(DFXLEN)- THE ARRAY SPECIFIES WHICH FIRST AND SECOND
C                   DERIVATIVES ARE TO BE COMPUTED BY CONTAINING THE
C                   NUMBER OF THE VARIABLES
C      IERR       - THE PARAMETER SHOWS THE REASON FOR TERMINATING THE
C                   SUBROUTINE. ON RETURN IERR COULD CONTAIN THE FOLLOW-
C                   ING VALUES:
C                   IERR = 0 : SUCCESSFUL TERMINATION.
C                   IERR > 0 : AN ERROR HAS BEEN DETECTED. FOR FURTHER
C                              INFORMATION CF. SUBROUTINE SYMERR.
C
C***********************************************************************
C
      INTEGER PIIS,PVIS,PIIC,PVIC,PIRC,PVRC,PIVA,PVVA,PIFN,PXFN,PVFN
C     INTEGER PVGR,PVHE,PVPF,PIV,LIIS,LVIS,LIIC,LVIC,LIRC,LVRC,LIVA
      INTEGER PVGR,PVHE,PVPF,LIIS,LVIS,LIIC,LVIC,LIRC,LVRC,LIVA
      INTEGER LVVA,LIFN,LVFN,LVGR,LVHE,LVPF,LIV,LGST,LHST
      INTEGER MPIIS,MPVIS,MPIIC,MPVIC,MPIRC,MPVRC,MPIVA,MPVVA,MPIFN
      INTEGER MPVFN,MPVPF,MPIV,MODE
C
      INTEGER IVAL
      DOUBLE PRECISION FVAL
      INTEGER I,J,K,L,PC
C
      INTEGER GSMDEP
      PARAMETER (GSMDEP=10)
C
      INTEGER GRAD
      PARAMETER (GRAD=1)
C
      IERR=0
C
      PIIS=IWA(1)
      PVIS=IWA(2)
      PIIC=IWA(3)
      PVIC=IWA(4)
      PIRC=IWA(5)
      PVRC=IWA(6)
      PIVA=IWA(7)
      PVVA=IWA(8)
      PIFN=IWA(9)
      PXFN=IWA(10)
      PVFN=IWA(11)
      PVGR=IWA(12)
      IF (PVGR .EQ. 0) THEN
        IERR=67
        RETURN
      ENDIF
      PVHE=1
      PVPF=IWA(14)
C     PIV=IWA(15)
      LIIS=16
      LVIS=LIIS+PIIS*5
      LIIC=LVIS+PVIS
      LVIC=LIIC+PIIC*4
      LIRC=LVIC+PVIC
      LVRC=1
      LIVA=LIRC+PIRC*4
      LVVA=LVRC+PVRC
      LIFN=LIVA+PIVA*3
      LVFN=LVVA+PVVA
      LVGR=LVFN+PVFN
      LVHE=LVGR+PVGR*PVVA
      LVPF=LIFN+PIFN*7
      LIV=LVPF+PVPF
      LGST=LVHE+PVHE*1
      LHST=LGST+GSMDEP*PVVA
      IF (N .NE. PVVA) THEN
        IERR=43
        RETURN
      ENDIF
      IF (M .NE. PVFN-(PIFN-PXFN)) THEN
        IERR=44
        RETURN
      ENDIF
      DO 10 I=1,N
        WA(LVVA+I-1)=X(I)
 10   CONTINUE
      K=0
      DO 20 I=1,PXFN
        IF (IWA(LIFN+(I-1)+(1-1)*PIFN) .EQ. 0) THEN
          K=K+1
          IF (ACTIVE(K)) THEN
            PC=IWA(LIFN+(I-1)+(6-1)*PIFN)
            MPIIS=MAX(1,IWA(1))
            MPVIS=MAX(1,IWA(2))
            MPIIC=MAX(1,IWA(3))
            MPVIC=MAX(1,IWA(4))
            MPIRC=MAX(1,IWA(5))
            MPVRC=MAX(1,IWA(6))
            MPIVA=MAX(1,IWA(7))
            MPVVA=MAX(1,IWA(8))
            MPIFN=MAX(1,IWA(9))
            MPVFN=MAX(1,IWA(11))
            MPVPF=MAX(1,IWA(14))
            MPIV =MAX(1,IWA(15))
            MODE=GRAD
            DO 30 L=LIV,LIV+MPIV-1
              IWA(L)=0
 30         CONTINUE
            CALL EVAL(PC,IVAL,FVAL,MODE,PVVA,1,
     1                MPIIS,MPVIS,MPIIC,MPVIC,MPIRC,MPVRC,MPIVA,MPVVA,
     2                MPIFN,MPVFN,PVHE,MPVPF,MPIV,IWA(LIIS),IWA(LVIS),
     3                IWA(LIIC),IWA(LVIC),IWA(LIRC),WA(LVRC),IWA(LIVA),
     4                WA(LVVA),IWA(LIFN),WA(LVFN),WA(LVGR),WA(LVHE),
     5                IWA(LVPF),IWA(LIV),WA(LGST),WA(LHST),DFX,DFXLEN,
     6                IERR)
            IF (IERR .NE. 0) RETURN
            F(K)=WA(LVFN+IWA(LIFN+(I-1)+(4-1)*PIFN)-1)
            DO 11 L=1,PVVA
              DF(K,L)=WA(LVGR+(IWA(LIFN+(I-1)+(5-1)*PIFN)-1)+(L-1)*PVGR)
 11         CONTINUE
          ENDIF
        ELSE IF (IWA(LIFN+(I-1)+(1-1)*PIFN) .EQ. 1) THEN
          DO 15 J=1,IWA(LIFN+(I-1)+(3-1)*PIFN)
            K=K+1
            IF (ACTIVE(K)) THEN
              PC=IWA(LIFN+(I-1)+(6-1)*PIFN)
              MPIIS=MAX(1,IWA(1))
              MPVIS=MAX(1,IWA(2))
              MPIIC=MAX(1,IWA(3))
              MPVIC=MAX(1,IWA(4))
              MPIRC=MAX(1,IWA(5))
              MPVRC=MAX(1,IWA(6))
              MPIVA=MAX(1,IWA(7))
              MPVVA=MAX(1,IWA(8))
              MPIFN=MAX(1,IWA(9))
              MPVFN=MAX(1,IWA(11))
              MPVPF=MAX(1,IWA(14))
              MPIV =MAX(1,IWA(15))
              MODE=GRAD
              DO 40 L=LIV,LIV+MPIV-1
                IWA(L)=0
 40           CONTINUE
              IWA(LIV+IWA(LIFN+(I-1)+(2-1)*PIFN)-1)=J
              CALL EVAL(PC,IVAL,FVAL,MODE,PVVA,1,
     1                MPIIS,MPVIS,MPIIC,MPVIC,MPIRC,MPVRC,MPIVA,MPVVA,
     2                MPIFN,MPVFN,PVHE,MPVPF,MPIV,IWA(LIIS),IWA(LVIS),
     3                IWA(LIIC),IWA(LVIC),IWA(LIRC),WA(LVRC),IWA(LIVA),
     4                WA(LVVA),IWA(LIFN),WA(LVFN),WA(LVGR),WA(LVHE),
     5                IWA(LVPF),IWA(LIV),WA(LGST),WA(LHST),DFX,DFXLEN,
     6                IERR)
              IF (IERR .NE. 0) RETURN
              F(K)=WA(LVFN+IWA(LIFN+(I-1)+(4-1)*PIFN)+(J-1)-1)
              DO 12 L=1,PVVA
                DF(K,L)=WA(LVGR+(IWA(LIFN+(I-1)+(5-1)*PIFN)-1+(J-1))+
     1              (L-1)*PVGR)
 12           CONTINUE
            ENDIF
 15       CONTINUE
        ENDIF
 20   CONTINUE
      RETURN
      END
C
C
C
      SUBROUTINE SYMHES (X,N,F,M,DF,DDF,MMAX,ACTIVE,WA,LWA,IWA,LIWA,
     1                   DFX,DFXLEN,IERR)
      INTEGER N,M,MMAX
      DOUBLE PRECISION X(N),F(M),DF(MMAX,N),DDF(MMAX,N*N)
      LOGICAL ACTIVE(M)
      INTEGER LWA,LIWA
      DOUBLE PRECISION WA(LWA)
      INTEGER IWA(LIWA)
      INTEGER DFXLEN
      INTEGER DFX(*)
      INTEGER IERR
C
C***********************************************************************
C
C   S Y M H E S   - EVALUATE SYMBOLICALLY DEFINED FUNCTIONS,
C                   CORRESPONDING GRADIENTS AND THE HESSIAN MATRIX
C
C   PARAMETERS:
C      X(N)          - ON INPUT, THE ONE-DIMENSIONAL ARRAY X HAS TO
C                      CONTAIN THE ARGUMENT THE FUNCTIONS ARE TO BE
C                      COMPUTED AT.
C      F(M)          - ON RETURN, F CONTAINS THE VALUES OF THE ACTIVE
C                      FUNCTIONS AT ARGUMENT X.
C      DF(MMAX,N)    - ON RETURN, DF CONTAINS THE GRADIENTS OF THE
C                      ACTIVE FUNCTIONS AT ARGUMENT X. IN THE DRIVING
C                      PROGRAM, THE ROW DIMENSION OF DF HAS TO BE EQUAL
C                      TO MMAX.
C      DDF(MMAX,N*N) - ON RETURN, DDF CONTAINS THE SECOND DERIVATIVES OF
C                      THE ACTIVE FUNCTIONS AT ARGUMENT X. IN THE
C                      DRIVING PROGRAM, THE ROW DIMENSION OF DDF HAS TO
C                      BE EQUAL TO MMAX
C      ACTIVE(M)     - THE LOGICAL ARRAY SPECIFIES WHICH OF THE M
C                      FUNCTIONS ARE TO BE COMPUTED (ACTIVE(K)=.TRUE.).
C      WA(LWA)       - REAL WORKING ARRAY, CONTAINS THE INTERMEDIATE
C                      CODE GENERATED BY SYMINP.
C      IWA(LIWA)     - INTEGER WORKING ARRAY, CF. WA(LWA).
C      DFX(DFXLEN)   - THE ARRAY SPECIFIES WHICH FIRST AND SECOND
C                      DERIVATIVES ARE TO BE COMPUTED BY CONTAINING THE
C                      NUMBER OF THE VARIABLES
C      IERR          - THE PARAMETER SHOWS THE REASON FOR TERMINATING
C                      THE SUBROUTINE. ON RETURN IERR COULD CONTAIN THE
C                      FOLLOWING VALUES:
C                      IERR = 0 : SUCCESSFUL TERMINATION.
C                      IERR > 0 : AN ERROR HAS BEEN DETECTED. FOR
C                                 FURTHER INFORMATION CF. SUBROUTINE
C                                 SYMERR.
C
C***********************************************************************
C
      INTEGER PIIS,PVIS,PIIC,PVIC,PIRC,PVRC,PIVA,PVVA,PIFN,PXFN,PVFN
C     INTEGER PVGR,PVHE,PVPF,PIV,LIIS,LVIS,LIIC,LVIC,LIRC,LVRC,LIVA
      INTEGER PVGR,PVHE,PVPF,LIIS,LVIS,LIIC,LVIC,LIRC,LVRC,LIVA
      INTEGER LVVA,LIFN,LVFN,LVGR,LVHE,LVPF,LIV,LGST,LHST
      INTEGER MPIIS,MPVIS,MPIIC,MPVIC,MPIRC,MPVRC,MPIVA,MPVVA,MPIFN
      INTEGER MPVFN,MPVPF,MPIV,MODE
C
      INTEGER IVAL
      DOUBLE PRECISION FVAL
      INTEGER I,J,K,L,P,PC
C
      INTEGER GSMDEP
      PARAMETER (GSMDEP=10)
C
      INTEGER HESS
      PARAMETER (HESS=2)
C
      IERR=0
C
      PIIS=IWA(1)
      PVIS=IWA(2)
      PIIC=IWA(3)
      PVIC=IWA(4)
      PIRC=IWA(5)
      PVRC=IWA(6)
      PIVA=IWA(7)
      PVVA=IWA(8)
      PIFN=IWA(9)
      PXFN=IWA(10)
      PVFN=IWA(11)
      PVGR=IWA(12)
      PVHE=IWA(13)
      IF (PVHE .EQ. 0) THEN
        IERR=62
        RETURN
      ENDIF
      PVPF=IWA(14)
C     PIV=IWA(15)
      LIIS=16
      LVIS=LIIS+PIIS*5
      LIIC=LVIS+PVIS
      LVIC=LIIC+PIIC*4
      LIRC=LVIC+PVIC
      LVRC=1
      LIVA=LIRC+PIRC*4
      LVVA=LVRC+PVRC
      LIFN=LIVA+PIVA*3
      LVFN=LVVA+PVVA
      LVGR=LVFN+PVFN
      LVHE=LVGR+PVGR*PVVA
      LVPF=LIFN+PIFN*7
      LIV=LVPF+PVPF
      LGST=LVHE+PVHE*PVVA*PVVA
      LHST=LGST+GSMDEP*PVVA
      IF (N .NE. PVVA) THEN
        IERR=43
        RETURN
      ENDIF
      IF (M .NE. PVFN-(PIFN-PXFN)) THEN
        IERR=44
        RETURN
      ENDIF
      DO 10 I=1,N
        WA(LVVA+I-1)=X(I)
 10   CONTINUE
      K=0
      DO 30 I=1,PXFN
        IF (IWA(LIFN+(I-1)+(1-1)*PIFN) .EQ. 0) THEN
          K=K+1
          IF (ACTIVE(K)) THEN
            PC=IWA(LIFN+(I-1)+(6-1)*PIFN)
            MPIIS=MAX(1,IWA(1))
            MPVIS=MAX(1,IWA(2))
            MPIIC=MAX(1,IWA(3))
            MPVIC=MAX(1,IWA(4))
            MPIRC=MAX(1,IWA(5))
            MPVRC=MAX(1,IWA(6))
            MPIVA=MAX(1,IWA(7))
            MPVVA=MAX(1,IWA(8))
            MPIFN=MAX(1,IWA(9))
            MPVFN=MAX(1,IWA(11))
            MPVPF=MAX(1,IWA(14))
            MPIV =MAX(1,IWA(15))
            MODE=HESS
            DO 40 L=LIV,LIV+MPIV-1
              IWA(L)=0
 40         CONTINUE
            CALL EVAL(PC,IVAL,FVAL,MODE,PVVA,PVVA,
     1                MPIIS,MPVIS,MPIIC,MPVIC,MPIRC,MPVRC,MPIVA,MPVVA,
     2                MPIFN,MPVFN,PVHE,MPVPF,MPIV,IWA(LIIS),IWA(LVIS),
     3                IWA(LIIC),IWA(LVIC),IWA(LIRC),WA(LVRC),IWA(LIVA),
     4                WA(LVVA),IWA(LIFN),WA(LVFN),WA(LVGR),WA(LVHE),
     5                IWA(LVPF),IWA(LIV),WA(LGST),WA(LHST),DFX,DFXLEN,
     6                IERR)
            IF (IERR .NE. 0) RETURN
            F(K)=WA(LVFN+IWA(LIFN+(I-1)+(4-1)*PIFN)-1)
            DO 11 L=1,PVVA
              DF(K,L)=WA(LVGR+(IWA(LIFN+(I-1)+(5-1)*PIFN)-1)+(L-1)*PVGR)
 11         CONTINUE
            DO 13 L=1,PVVA
              DO 12 P=1,PVVA
                DDF(K,(L-1)*PVVA+P)=WA(LVHE+(IWA(LIFN+(I-1)+(7-1)*PIFN)
     1              -1)+((L-1)*PVVA+P-1)*PVHE)
 12           CONTINUE
 13         CONTINUE
          ENDIF
        ELSE IF (IWA(LIFN+(I-1)+(1-1)*PIFN) .EQ. 1) THEN
          DO 20 J=1,IWA(LIFN+(I-1)+(3-1)*PIFN)
            K=K+1
            IF (ACTIVE(K)) THEN
              PC=IWA(LIFN+(I-1)+(6-1)*PIFN)
              MPIIS=MAX(1,IWA(1))
              MPVIS=MAX(1,IWA(2))
              MPIIC=MAX(1,IWA(3))
              MPVIC=MAX(1,IWA(4))
              MPIRC=MAX(1,IWA(5))
              MPVRC=MAX(1,IWA(6))
              MPIVA=MAX(1,IWA(7))
              MPVVA=MAX(1,IWA(8))
              MPIFN=MAX(1,IWA(9))
              MPVFN=MAX(1,IWA(11))
              MPVPF=MAX(1,IWA(14))
              MPIV =MAX(1,IWA(15))
              MODE=HESS
              DO 50 L=LIV,LIV+MPIV-1
                IWA(L)=0
 50           CONTINUE
              IWA(LIV+IWA(LIFN+(I-1)+(2-1)*PIFN)-1)=J
              CALL EVAL(PC,IVAL,FVAL,MODE,PVVA,PVVA,
     1                MPIIS,MPVIS,MPIIC,MPVIC,MPIRC,MPVRC,MPIVA,MPVVA,
     2                MPIFN,MPVFN,PVHE,MPVPF,MPIV,IWA(LIIS),IWA(LVIS),
     3                IWA(LIIC),IWA(LVIC),IWA(LIRC),WA(LVRC),IWA(LIVA),
     4                WA(LVVA),IWA(LIFN),WA(LVFN),WA(LVGR),WA(LVHE),
     5                IWA(LVPF),IWA(LIV),WA(LGST),WA(LHST),DFX,DFXLEN,
     6                IERR)
              IF (IERR .NE. 0) RETURN
              F(K)=WA(LVFN+IWA(LIFN+(I-1)+(4-1)*PIFN)+(J-1)-1)
              DO 14 L=1,PVVA
                DF(K,L)=WA(LVGR+(IWA(LIFN+(I-1)+(5-1)*PIFN)-1+(J-1))+
     1              (L-1)*PVGR)
 14           CONTINUE
              DO 16 L=1,PVVA
                DO 15 P=1,PVVA
                  DDF(K,(L-1)*PVVA+P)=WA(LVHE+(IWA(LIFN+(I-1)+
     1                  (7-1)*PIFN)-1+(J-1))+((L-1)*PVVA+P-1)*PVHE)
 15             CONTINUE
 16           CONTINUE
            ENDIF
 20       CONTINUE
        ENDIF
 30   CONTINUE
      RETURN
      END

