C ** Correction report.
C ** Correction 1. 05/02/93: 2 lines added ** 
C ** Correction 2. 18/07/94: 1 line corrected ** 
C ** Correction 3. 18/07/94: 1 line corrected ** 
C ** End of Correction report.
C  THIS VERSION: 18/07/1994 AT 09:23:34 AM.
C
C  ** FOR THE CRAY 2, LINES STARTING 'CDIR$ IVDEP' TELL THE COMPILER TO
C     IGNORE POTENTIAL VECTOR DEPENDENCIES AS THEY ARE KNOWN TO BE O.K.
C
CS    SUBROUTINE SINITW( N , NG, NEL   , IELING, LELING, ISTADG, LSTADG,
CD    SUBROUTINE DINITW( N , NG, NEL   , IELING, LELING, ISTADG, LSTADG,
     *                   IELVAR, LELVAR, ISTAEV, LSTAEV, INTVAR, LNTVAR,
     *                   ISTADH, LSTADH, ICNA  , LICNA , ISTADA, LSTADA,
     *                   GXEQX , LGXEQX, INTREP, LINTRE, LFUVAL, ALTRIV, 
     *                   DIRECT, FDGRAD, LFXI  , LGXI  , LHXI  , LGGFX , 
     *                   LDX   , LGRJAC, LQGRAD, LBREAK, LP    , LXCP  , 
     *                   LX0   , LGX0  , LDELTX, LBND  , LWKSTR, LSPTRS, 
     *                   LSELTS, LINDEX, LSWKSP, LSTAGV, LSTAJC, LIUSED, 
     *                   LFREEC, LNNONZ, LNONZ2, LSYMMD, LSYMMH, LSLGRP, 
     *                   LSVGRP, LGCOLJ, LVALJR, LSEND , LNPTRS, LNELTS, 
     *                   LNNDEX, LNWKSP, LNSTGV, LNSTJC, LNIUSE, LNFREC, 
     *                   LNNNON, LNNNO2, LNSYMD, LNSYMH, LNLGRP, LNVGRP, 
     *                   LNGCLJ, LNVLJR, LNQGRD, LNBRAK, LNP   , LNBND , 
     *                   LNFXI , LNGXI , LNGUVL, LNHXI , LNHUVL, LNGGFX, 
     *                   LNDX  , LNGRJC, LIWK2 , LWK2  , MAXSIN, NINVAR, 
     *                   NTYPE , NSETS , MAXSEL, LSTYPE, LSSWTR, LSSIWT, 
     *                   LSIWTR, LSWTRA, LNTYPE, LNSWTR, LNSIWT, LNIWTR, 
     *                   LNWTRA, LSISET, LSSVSE, LNISET, LNSVSE, RANGES,
     *                   IWK   , LIWK  , WK    , LWK   ,
     *                   IPRINT, IOUT  , INFORM )
C
C  --------------------------------------------------------------------
C
C  COMPUTE THE STARTING ADDRESSES FOR THE PARTITIONS OF THE WORKSPACE
C  ARRAYS FUVALS, IWK AND WK. ALSO FILL RELEVANT PORTIONS OF IWK.
C  THE ADDRESSES FOR IWK AND WK ARE DESCRIBED HERE; THOSE FOR FUVALS
C  ARE DESCRIBED IN THE INTRODUCTORY COMMENTS TO SUBROUTINE SBMIN.
C
C  NICK GOULD, 20TH JUNE 1990.
C  FOR CGT PRODUCTIONS.
C
C  --------------------------------------------------------------------
C
      INTEGER          N , NG, NEL   , INFORM, LWK   , LIWK
      INTEGER          LFUVAL, LELVAR, LSTAEV, LSTADH, LELING, LNTVAR
      INTEGER          LGXEQX, LSTADG, LICNA , LSTADA, IOUT  , IPRINT
      INTEGER          LFXI  , LGXI  , LHXI  , LGGFX , LDX   , LGRJAC
      INTEGER          LQGRAD, LBREAK, LP    , LXCP  , LBND  , LWKSTR
      INTEGER          LGX0  , LDELTX, LSPTRS, LSELTS, LINDEX, LX0
      INTEGER          LSWKSP, LSTAGV, LSTAJC, LIUSED, LFREEC, LINTRE
      INTEGER          LNNONZ, LNONZ2, LSYMMD, LSYMMH, LSISET, LSSVSE
      INTEGER          LSLGRP, LSVGRP, LGCOLJ, LVALJR, LSEND
      INTEGER          LNPTRS, LNELTS, LNNDEX, LNWKSP, LNSTGV, LNSTJC
      INTEGER          LNIUSE, LNNNON, LNNNO2, LNSYMD, LNSYMH, LNLGRP
      INTEGER          LNVGRP, LNGCLJ, LNVLJR, LNFREC, LNISET, LNSVSE
      INTEGER          LNQGRD, LNBRAK, LNP   , LNBND,  LNFXI , LNGXI
      INTEGER          LNGUVL, LNHXI , LNHUVL, LNGGFX, LNDX  , LNGRJC
      INTEGER          LIWK2 , LWK2  , NINVAR, MAXSEL, MAXSIN, NTYPE 
      INTEGER          LSTYPE, LSSWTR, LSSIWT, LSIWTR, LSWTRA, NSETS
      INTEGER          LNTYPE, LNSWTR, LNSIWT, LNIWTR, LNWTRA
      EXTERNAL         RANGES 
      LOGICAL          DIRECT, ALTRIV, FDGRAD
      INTEGER          IELVAR( LELVAR       ), ISTAEV( LSTAEV       )
      INTEGER          ISTADH( LSTADH       ), IWK   ( LIWK         )
      INTEGER          INTVAR( LNTVAR       ), ISTADG( LSTADG       )
      INTEGER          ICNA  ( LICNA        ), ISTADA( LSTADA       )
      INTEGER          IELING( LELING       )
CS    REAL             WK    ( LWK          )
CD    DOUBLE PRECISION WK    ( LWK          )
      LOGICAL          GXEQX ( LGXEQX       ), INTREP( LINTRE       )
C
C  LOCAL VARIABLES.
C
      INTEGER          I , J , K , L , NVARGP, IIELTS, IENTRY, IG, IS
      INTEGER          LEND  , NSIZEH, NEL1  , NGEL  , NG1   , IEL
      INTEGER          LWFREE, LIWFRE, NELVR , LW1   , LIWFRO, LWFREO
      INTEGER          IELL  , ITYPE , ISOFAR, ISTARJ, IVARP1, IVAR
      INTEGER          JSET  , INEXT , NEWVAR, NEWSET, IPT   , ISTRT
      INTEGER          NINVR , II, JJ, KK, LL
      LOGICAL          NONTRV, ALLLIN, VRUSED          
CS    REAL             ZERO  , ONE   
CD    DOUBLE PRECISION ZERO  , ONE   
CS    PARAMETER      ( ZERO = 0.0E+0 , ONE = 1.0E+0  )
CD    PARAMETER      ( ZERO = 0.0D+0 , ONE = 1.0D+0  )
C
C  EXTERNAL SUBROUTINES AND FUNCTIONS USED.
C
CS    EXTERNAL         SSYMMH
CD    EXTERNAL         DSYMMH
C
C  INTRINSIC FUNCTIONS.
C
      INTRINSIC        MAX
C
C  SET CONSTANTS.
C
      NEL1   = NEL + 1
      NG1    = NG  + 1
      NGEL   = ISTADG( NG1 ) - 1
      ALLLIN = NEL .EQ. 0
C
C  SET UP INTVAR, THE STARTING ADDRESSES FOR THE ELEMENT GRADIENTS
C  WITH RESPECT TO THEIR INTERNAL VARIABLES. ALSO COMPUTE MAXSIN,
C  THE MAXIMUM NUMBER OF INTERNAL VARIABLES IN AN ELEMENT.
C
      IF ( .NOT. ALLLIN ) THEN
         K                = INTVAR( 1 )
         MAXSIN           = K
         INTVAR( 1 )      = NEL1
         DO 10 IEL        = 2, NEL
            L             = INTVAR( IEL )
            INTVAR( IEL ) = INTVAR( IEL - 1 ) + K
            K             = L
            MAXSIN        = MAX( MAXSIN, K )
   10    CONTINUE
         INTVAR( NEL1 ) = INTVAR( NEL ) + K
      ELSE
         INTVAR( 1 ) = 1
         MAXSIN      = 0
      END IF
C
C  COMPUTE THE TOTAL NUMBER OF INTERNAL VARIABLES.
C
      NINVAR = INTVAR( NEL1 ) - INTVAR( 1 )
C
C  CALCULATE THE LENGTH, IIELTS, OF WORKSPACE REQUIRED TO
C  DETERMINE WHICH ELEMENTS USE EACH OF THE VARIABLES.
C  ALSO FIND THE MAXIMUM NUMBER OF VARIABLES IN AN ELEMENT, MAXSEL.
C  THIS IS A DUMMY RUN FOR LOOP 130 MERELY TO CALCULATE THE 
C  SPACE REQUIRED.
C
      IF ( LIWK .LT. N ) THEN
         WRITE( IOUT, 2030 ) N - LIWK
         INFORM = 4
         RETURN
      END IF
C
C  IWK( I ) WILL BE USED AS A LIST OF LINKS CHAINING THE ELEMENTS USING
C  VARIABLE I. IF IWK( I ) IS NEGATIVE, THE LIST IS EMPTY.
C
CDIR$ IVDEP
      DO 20 I     = 1, N
         IWK( I ) = - 1
   20 CONTINUE   
      IIELTS = N
      MAXSEL = 0
      IF ( .NOT. ALLLIN ) THEN
C
C  LOOP OVER THE GROUPS, CONSIDERING EACH NONLINEAR ELEMENT IN TURN.
C
         DO 50 I   = 1, NGEL
            IEL    = IELING( I )
            MAXSEL = MAX( MAXSEL, ISTAEV( IEL + 1 ) - ISTAEV( IEL ) )
C
C  LOOP ON THE VARIABLES FROM THE I-TH ELEMENT.
C
            DO 40 K   = ISTAEV( IEL ), ISTAEV( IEL + 1 ) - 1
               IENTRY = IELVAR( K )
               IF ( IWK( IENTRY ) .GE. 0 ) THEN
C
C  IF WE HAVE REACHED THE END OF THE LIST OF THE ELEMENTS USING 
C  THE VARIABLE IELVAR( K ), ADD THE IEL-TH ELEMENT TO IT.
C  OTHERWISE, FIND THE NEXT ENTRY IN THE LIST.
C
   30             CONTINUE
                  IF ( IWK( IENTRY ) .GT. 0 ) THEN
                     IENTRY = IWK( IENTRY )
                     GO TO 30
                  ELSE
                     IIELTS = IIELTS + 1
                     IF ( IIELTS .GT. LIWK ) THEN
                        WRITE( IOUT, 2030 ) IIELTS - LIWK
                        INFORM = 4
                        RETURN
                     END IF
                     IWK( IENTRY ) = IIELTS
                     IWK( IIELTS ) = 0
                  END IF
               ELSE
C
C  THE LIST OF ELEMENTS INVOLVING THE VARIABLE IELVAR( K ) WAS 
C  PREVIOUSLY EMPTY. INDICATE THAT THE LIST HAS NOW BEEN STARTED AND 
C  THAT ITS END HAS BEEN REACHED.
C
                  IWK( IENTRY ) = 0
               END IF
   40       CONTINUE
   50    CONTINUE
      END IF
C
C  -- CALCULATE THE STARTING ADDRESSES FOR THE INTEGER WORKSPACE. --
C
C  IWK( LSPTRS + J ), J = 1, ..., IIELTS,  WILL CONTAIN THE LINKS FOR
C  THE LISTS OF NONLINEAR ELEMENTS WHICH USE EACH VARIABLE.
C
      LSPTRS = 0
C
C  IWK( LSELTS + J ), J = 1, ..., IIELTS, WILL CONTAIN THE LISTS OF 
C  NONLINEAR ELEMENTS CORRESPONDING TO THE PREVIOUSLY MENTIONED LINKS.
C
      LSELTS = LSPTRS + IIELTS
C
C  IWK( LINDEX + J ), J = 1, ..., N, WILL CONTAIN THE STATUS OF THE
C  J-TH VARIABLE AS THE CURRENT ITERATION PROGRESSES. POSSIBLE VALUES
C  ARE 0 IF THE VARIABLE LIES AWAY FROM ITS BOUNDS, 1 AND 2 IF IT LIES
C  ON ITS LOWER OR UPPER BOUNDS (RESPECTIVELY) - THESE MAY BE PROBLEM
C  BOUNDS OR TRUST REGION BOUNDS, AND 3 IF THE VARIABLE IS FIXED.
C
      LINDEX = LSELTS + IIELTS
C
C  IWK( LSWKSP + J ), J = 1, ..., MAX( NGEL, NEL, N + N ), IS USED FOR
C  WORKSPACE BY THE MATRIX-VECTOR PRODUCT SUBROUTINE HESPRD.
C
      LSWKSP = LINDEX + N
C
C  IWK( LIUSED + J ), J = 1, ..., MAX( N, NG ) WILL BE USED AS 
C  WORKSPACE BY THE MATRIX-VECTOR PRODUCT SUBROUTINE HESPRD.
C
      IF ( DIRECT ) THEN
         LIUSED = LSWKSP + MAX( NGEL, NEL, N + N )
      ELSE
         LIUSED = LSWKSP + MAX( NGEL, NEL, N )
      END IF
C
C  IWK( LFREED + J ), J = 1, ..., NFREEC WILL GIVE THE INDICES OF THE
C  VARIABLES WHICH ARE CONSIDERED TO BE FREE FROM THEIR BOUNDS AT THE
C  CURRENT GENERALIZED CAUCHY POINT.
C
      LFREEC = LIUSED + MAX( N, NG )
C
C  IWK( LNNONZ + J ), J = 1, ..., NNNONZ WILL GIVE THE INDICES OF THE
C  NONZEROS IN THE VECTOR OBTAINED AS A RESULT OF THE MATRIX-VECTOR
C  PRODUCT FROM SUBROUTINE HESPRD.
C
      LNNONZ = LFREEC + N
C
C  IWK( LNONZ2 + J ), J = 1, ..., NG, WILL BE USED AS FURTHER
C  WORKSPACE BY THE MATRIX-VECTOR PRODUCT SUBROUTINE HESPRD.
C
      LNONZ2 = LNNONZ + N
C
C  IWK( LSYMMD + J ), J = 1, ..., MAXSIN, WILL GIVE THE LOCATION OF
C  THE J-TH DIAGONAL OF A MAXSIN BY MAXSIN SYMMETRIC MATRIX IN AN 
C  UPPER TRIANGULAR STORAGE SCHEME.
C
      LSYMMD = LNONZ2 + NG
C
C  IWK( LSYMMH + ( I - 1 ) * MAXSIN + J ), I, J = 1, ..., MAXSIN, 
C  WILL GIVE THE LOCATION OF THE (I,J)-TH ENTRY OF A MAXSIN BY 
C  MAXSIN SYMMETRIC MATRIX IN AN UPPER TRIANGULAR STORAGE SCHEME.
C
      LSYMMH = LSYMMD + MAXSIN
C
C  IWK( LSLGRP + J ), J = 1, ..., NGEL, WILL CONTAIN THE NUMBER OF 
C  THE GROUP WHICH USES NONLINEAR ELEMENT J.
C
      LSLGRP = LSYMMH + MAXSIN * MAXSIN
C
C  IWK( LSTAJC + J ), J = 1, ..., N, WILL CONTAIN THE STARTING 
C  ADDRESSES FOR THE LIST OF NONTRIVIAL GROUPS WHICH USE THE 
C  J-TH VARIABLE. IWK( LSTAJC + N + 1 ) WILL POINT TO THE FIRST FREE 
C  LOCATION IN IWK AFTER THE LIST OF NONTRIVIAL GROUPS FOR THE 
C  N-TH VARIABLE.
C
      LSTAJC = LSLGRP + NGEL
C
C  IWK( LSTAGV + J ), J = 1, ..., NG, WILL CONTAIN THE STARTING 
C  ADDRESSES FOR THE LIST OF VARIABLES WHICH OCCUR IN THE J-TH GROUP. 
C  IWK( LSTAGV + NG + 1 ) WILL POINT TO THE FIRST FREE LOCATION
C  IN IWK AFTER THE LIST OF VARIABLES FOR THE NG-TH GROUP.
C
      LSTAGV = LSTAJC + N + 1
C
C  IWK( LSVGRP + J ), J = 1, ..., NVARGP, WILL CONTAIN THE INDICES OF
C  THE VARIABLES WHICH ARE USED BY EACH GROUP IN TURN. THOSE FOR GROUP I
C  OCCUR IN LOCATIONS IWK( LSTAGV + I ) TO IWK( LSTAGV + I + 1 ) - 1.
C
      LSVGRP = LSTAGV + NG + 1
C
C  CHECK THAT THERE IS SUFFICIENT WORKSPACE.
C
      IF ( LSVGRP .GT. LIWK ) THEN
         WRITE( IOUT, 2030 ) LSVGRP - LIWK
         INFORM = 4
         RETURN
      END IF
C
C  IWK( LGCOLJ + J ), J = 1, ..., NVARGP, WILL CONTAIN THE INDICES OF
C  THE NONTRIVIAL GROUPS WHICH USE EACH VARIABLE IN TURN. THOSE FOR 
C  VARIABLE I OCCUR IN LOCATIONS IWK( LSTAJC + I ) TO 
C  IWK( LSTAJC + I + 1 ) - 1.
C
      LGCOLJ = LSVGRP + 1
C
C  DETERMINE WHICH ELEMENTS USE EACH VARIABLE. INITIALIZATION.
C
      IF ( .NOT. ALLLIN ) THEN
C
C  IWK( I ) WILL BE USED AS A LIST OF LINKS CHAINING THE ELEMENTS USING
C  VARIABLE I. IF IWK( I ) IS NEGATIVE, THE LIST IS EMPTY.
C
CDIR$ IVDEP
         DO 100 I      = 1, N
            IWK( I ) = - 1
  100    CONTINUE
         IIELTS   = N
C
C  LOOP OVER THE GROUPS, CONSIDERING EACH NONLINEAR ELEMENT IN TURN.
C
         DO 130 I = 1, NGEL
            IEL   = IELING( I )
C
C  LOOP ON THE VARIABLES OF THE I-TH ELEMENT.
C
            DO 120 K  = ISTAEV( IEL ), ISTAEV( IEL + 1 ) - 1
               IENTRY = IELVAR( K )
               IF ( IWK( IENTRY ) .GE. 0 ) THEN
C
C  IF WE HAVE REACHED THE END OF THE LIST OF THE ELEMENTS USING 
C  THE VARIABLE IELVAR( K ), ADD THE I-TH ELEMENT TO IT AND 
C  RECORD THAT THE END OF THE LIST HAS OCCURED.
C  OTHERWISE, FIND THE NEXT ENTRY IN THE LIST.
C
  110             CONTINUE
                  IF ( IWK( IENTRY ) .GT. 0 ) THEN
                     IENTRY = IWK( IENTRY )
                     GO TO 110
                  ELSE
                     IIELTS                 = IIELTS + 1
                     IWK( IENTRY )          = IIELTS
                     IWK( LSELTS + IIELTS ) = I
C                    IWK( LSELTS + IIELTS ) = IEL
                     IWK( IIELTS )          = 0
                  END IF
               ELSE
C
C  THE LIST OF ELEMENTS INVOLVING THE VARIABLE IELVAR( K ) WAS 
C  PREVIOUSLY EMPTY. INDICATE THAT THE LIST HAS NOW BEEN STARTED, 
C  RECORD THE ELEMENT WHICH CONTAINS THE VARIABLE AND INDICATE
C  THAT THE END OF THE LIST HAS BEEN REACHED.
C
                  IWK( LSELTS + IENTRY ) = I
C                 IWK( LSELTS + IENTRY ) = IEL
                  IWK( IENTRY )          = 0
               END IF
  120       CONTINUE
  130    CONTINUE
      END IF
C
C  SET UP SYMMETRIC ADDRESSES FOR THE UPPER TRIANGULAR STORAGE
C  SCHEMES FOR THE ELEMENT HESSIANS.
C
CS    IF ( MAXSIN .GT. 0 ) CALL SSYMMH ( MAXSIN, IWK( LSYMMH + 1 ),
CD    IF ( MAXSIN .GT. 0 ) CALL DSYMMH ( MAXSIN, IWK( LSYMMH + 1 ),
     *                                  IWK( LSYMMD + 1 ) )
C
C  SET UP THE STARTING ADDRESSES FOR THE ELEMENT HESSIANS
C  WITH RESPECT TO THEIR INTERNAL VARIABLES AND A POINTER BEYOND
C  THE END OF THE SPACE REQUIRED FOR THE HESSIANS. 
C
      LGGFX = INTVAR( NEL1 )
      IF ( .NOT. ALLLIN ) THEN
         DO 140 I       = 1, NEL
            ISTADH( I ) = LGGFX
            NSIZEH      = INTVAR( I + 1 ) - INTVAR( I )
            LGGFX       = LGGFX + NSIZEH * ( NSIZEH + 1 ) / 2
  140    CONTINUE
      END IF
      ISTADH( NEL1 ) = LGGFX
C
C  -- CALCULATE THE STARTING ADDRESSES FOR THE REAL WORKSPACE. --
C
C  WK( LQGRAD + J ), J = 1, ..., N, WILL CONTAIN THE GRADIENT OF
C  THE QUADRATIC MODEL AT THE CURRENT ESTIMATE OF THE MINIMIZER.
C
      LQGRAD = MAX( MAX( NG, MAXSEL ) + 2 * MAXSIN, 
     *              N + NINVAR + MAX( MAXSEL, NINVAR ) )
C
C  WK( LBREAK + J ), J = 1, ..., N, WILL CONTAIN THE BREAKPOINTS
C  ALONG THE CAUCHY ARC FROM THE CURRENT ESTIMATE OF THE MINIMIZER.
C
      LBREAK = LQGRAD + N
C
C  WK( LP + J ), J = 1, ..., N, WILL CONTAIN THE VECTOR REQUIRED
C  BY THE MATRIX-VECTOR PRODUCT SUBROUTINE HESPRD.
C
      LP = LBREAK + N
C
C  WK( LXCP + J ), J = 1, ..., N, WILL CONTAIN THE CURRENT 
C  (APPROXIMATE) CAUCHY POINT.
C
      LXCP = LP + N
C
C  WK( LX0 + J ), J = 1, ..., N, WILL CONTAIN THE CURRENT 
C  START OF THE CAUCHY SEARCH. THIS FEATURE IS ONLY USED IF
C  MORE THAN ONE CYCLE IS USED TO SOLVE THE BQP ACCURATELY.
C
      LX0 = LXCP + N
C
C  WK( LGX0 + J ), J = 1, ..., N, WILL CONTAIN THE GRADIENT AT THE 
C  CURRENT START OF THE CAUCHY SEARCH. THIS FEATURE IS ONLY USED IF
C  MORE THAN ONE CYCLE IS USED TO SOLVE THE BQP ACCURATELY.
C
      LGX0 = LX0 + N
C
C  WK( LDELTX + J ), J = 1, ..., N, WILL CONTAIN THE STEP TAKEN 
C  DURING A SINGLE CYCLE, IF MORE THAN ONE CYCLE IS USED TO SOLVE
C  THE BQP ACCURATELY.
C
      LDELTX = LGX0 + N
C
C  WK( LBND + 2 * ( J - 1 ) + I ), J = 1, ..., N, I = 1, 2, WILL
C  CONTAIN THE CURRENT LOWER (I=1) AND UPPER (I=2) BOUNDS ON THE 
C  VARIABLES DEFINED BY THE INTERSECTION OF THE TRUST REGION WITH 
C  THE FEASIBLE BOX.
C
      LBND = LDELTX + N
C
C  WK( LWKSTR + J ), J = 1, ..., LWK2, IS THE REMAINING REAL 
C  WORKSPACE WHICH IS FREE FOR OTHER PURPOSES, SUCH AS FOR FORMING
C  THE FACTORIZATION OF THE MODEL HESSIAN, IF REQUIRED.
C
      LWKSTR = LBND + N + N
      LWK2   = LWK  - LWKSTR
C
C  CHECK THAT THERE IS SUFFICIENT REAL WORKSPACE.
C
      IF ( LWKSTR .GT. LWK ) THEN
         INFORM = 5
         WRITE( IOUT, 2040 ) LWKSTR - LWK
         RETURN
      END IF
C
C  SET THE LENGTH OF EACH PARTITION OF THE REAL WORKSPACE ARRAY 
C  FUVALS FOR ARRAY BOUND CHECKING IN CALLS TO OTHER SUBPROGRAMS.
C
      LNQGRD = LBREAK - LQGRAD
      LNBRAK = LP     - LBREAK
      LNP    = LXCP   - LP
      LNBND  = LWKSTR - LBND
C
C  STORE THE INDICES OF VARIABLES WHICH APPEARS IN EACH GROUP
C  AND HOW MANY GROUPS USE EACH VARIABLE. START BY INITIALIZING
C  COUNTING ARRAYS TO ZERO.
C
CDIR$ IVDEP
      DO 150 J                 = 1, N
         IWK( LSWKSP + J )     = 0
         IWK( LSTAJC + J + 1 ) = 0
  150 CONTINUE
C
C  ALTRIV SPECIFIES WHETHER ALL THE GROUPS ARE TRIVIAL.
C 
      ALTRIV = .TRUE.
C
C  COUNT THE TOTAL NUMBER OF VARIABLES IN ALL THE GROUPS, NVARGP.
C 
      NVARGP            = 0
      IWK( LSTAGV + 1 ) = 1
C
C  LOOP OVER THE GROUPS. SEE IF THE IG-TH GROUP IS TRIVIAL. 
C
      DO 200 IG  = 1, NG
         NONTRV  = .NOT. GXEQX( IG )
C
C  CHECK TO SEE IF ALL OF THE GROUPS ARE TRIVIAL.
C
         IF ( NONTRV ) ALTRIV = .FALSE.
C
C  LOOP OVER THE NONLINEAR ELEMENTS FROM THE IG-TH GROUP.
C
         DO 170 K = ISTADG( IG ), ISTADG( IG + 1 ) - 1
            IEL   = IELING( K )
C
C  RUN THROUGH ALL THE ELEMENTAL VARIABLES CHANGING THE I-TH ENTRY OF
C  IWK( LSWKSP ) FROM ZERO TO ONE IF VARIABLE I APPEARS IN AN ELEMENT.
C
            DO 160 J = ISTAEV( IEL ), ISTAEV( IEL + 1 ) - 1
               I     = IELVAR( J )
               IF ( IWK( LSWKSP + I ) .EQ. 0 ) THEN
                  IWK( LSWKSP + I ) = 1
C
C  IF THERE IS SUFFICIENT ROOM, RECORD THE NONLINEAR VARIABLES FROM 
C  THE IG-TH GROUP.
C
                  IF ( LGCOLJ .GT. LIWK ) THEN
                     WRITE( IOUT, 2030 ) LGCOLJ - LIWK
                     INFORM = 4
                     RETURN
                  END IF
                  IWK( LGCOLJ ) = I
                  LGCOLJ        = LGCOLJ + 1
                  NVARGP        = NVARGP + 1
               END IF
  160       CONTINUE
C
C  RECORD THAT NONLINEAR ELEMENT K OCCURS IN GROUP IELGRP(IEL).
C
            IWK( LSLGRP + K ) = IG
  170    CONTINUE
C
C  CONSIDER VARIABLES WHICH ARISE FROM THE LINEAR ELEMENT.
C
         DO 180 J = ISTADA( IG ), ISTADA( IG + 1 ) - 1
            I     = ICNA( J )
            IF ( IWK( LSWKSP + I ) .EQ. 0 ) THEN
               IWK( LSWKSP + I ) = 1
C
C  IF THERE IS SUFFICIENT ROOM, RECORD THE LINEAR VARIABLES FROM 
C  THE IG-TH GROUP.
C
               IF ( LGCOLJ .GT. LIWK ) THEN
                  WRITE( IOUT, 2030 ) LGCOLJ - LIWK
                  INFORM = 4
                  RETURN
               END IF
               IWK( LGCOLJ ) = I
               LGCOLJ        = LGCOLJ + 1
               NVARGP        = NVARGP + 1
            END IF
  180    CONTINUE
C
C  RESET THE STATUS ARRAY IWK( LSWKSP ) TO ZERO.
C
         DO 190 J = IWK( LSTAGV + IG ), NVARGP
            L     = IWK( LSVGRP + J )
            IWK( LSWKSP + L ) = 0
C
C  RECORD THAT ONE FURTHER NONTRIVIAL GROUP USES VARIABLE L.
C
            IF ( NONTRV ) IWK( LSTAJC + L + 1 ) = 
     *                    IWK( LSTAJC + L + 1 ) + 1               
  190    CONTINUE
C
C  RECORD THE STARTING ADDRESS OF THE VARIABLES IN THE NEXT GROUP.
C  
         IWK( LSTAGV + IG + 1 ) = NVARGP + 1
  200 CONTINUE
C
C  -- CONTINUE SETTING STARTING ADDRESSES FOR PARTITIONS OF IWK. --
C
C  IWK( LVALJR + J ), J = 1, ..., NVARGP, WILL CONTAIN THE POSITIONS IN
C  FUVALS (RELATIVE TO THE STARTING ADDRESS LGRJAC) OF THE NONZEROS OF
C  THE JACOBIAN OF THE GROUPS CORRESPONDING TO THE VARIABLES AS ORDERED
C  IN IWK( LSVGRP + J ).
C
      LVALJR = LGCOLJ + NVARGP
C
C  LSEND GIVES THE TOTAL FIXED AMOUNT OF INTEGER WORKSPACE USED. LIWK2
C  GIVES THE AMOUNT OF WORKSPACE WHICH IS FREE FOR OTHER PURPOSES, SUCH
C  AS FOR FORMING THE FACTORIZATION OF THE MODEL HESSIAN, IF REQUIRED.
C  
      LSEND  = LVALJR + NVARGP
C
C  CHECK THAT THERE IS SUFFICIENT INTEGER WORKSPACE.
C
      IF ( LSEND .GT. LIWK ) THEN
         WRITE( IOUT, 2030 ) LSEND - LIWK
         INFORM = 4
         RETURN
      END IF
C
C  SET THE LENGTH OF EACH PARTITION OF THE INTEGER WORKSPACE FOR
C  ARRAY BOUND CHECKING IN CALLS TO OTHER SUBPROGRAMS.
C
      LNPTRS = MAX( 1, LSELTS - LSPTRS )
      LNELTS = MAX( 1, LINDEX - LSELTS )
      LNNDEX = MAX( 1, LSWKSP - LINDEX )
      LNWKSP = MAX( 1, LIUSED - LSWKSP )
      LNIUSE = MAX( 1, LFREEC - LIUSED )
      LNFREC = MAX( 1, LNNONZ - LFREEC )
      LNNNON = MAX( 1, LNONZ2 - LNNONZ )
      LNNNO2 = MAX( 1, LSYMMD - LNONZ2 )
      LNSYMD = MAX( 1, LSYMMH - LSYMMD )
      LNSYMH = MAX( 1, LSLGRP - LSYMMH )
      LNLGRP = MAX( 1, LSTAJC - LSLGRP )
      LNSTJC = MAX( 1, LSTAGV - LSTAJC )
      LNSTGV = MAX( 1, LSVGRP - LSTAGV )
      LNVGRP = MAX( 1, LGCOLJ - LSVGRP )
      LNGCLJ = MAX( 1, LVALJR - LGCOLJ )
      LNVLJR = MAX( 1, LSEND  - LVALJR )
C
C  SET THE STARTING ADDRESSES FOR THE LISTS OF NONTRIVIAL GROUPS 
C  WHICH USE EACH VARIABLE IN TURN.
C
      K           = LSTAJC + 1
      IWK( K )    = 1
      DO 210 I    = 2, N + 1
         K        = K + 1
         IWK( K ) = IWK( K ) + IWK( K - 1 )
  210 CONTINUE
C
C  CONSIDER THE IG-TH GROUP IN ORDER TO ASSOCIATE VARIABLES WITH GROUPS.
C
      DO 230 IG = 1, NG
         IF ( .NOT. GXEQX( IG ) ) THEN
            DO 220 I = IWK( LSTAGV + IG ), IWK( LSTAGV + IG + 1 ) - 1
                L   = LSTAJC + IWK( LSVGRP + I )
C
C  RECORD THAT GROUP IG USES VARIABLE IWK( LSVGRP + I ).
C
               J                 = IWK( L )
               IWK( LGCOLJ + J ) = IG
C
C  STORE THE LOCATIONS IN THE JACOBIAN OF THE GROUPS OF THE NONZEROS
C  CORRESPONDING TO EACH VARIABLE IN THE IG-TH GROUP. INCREMENT THE
C  STARTING ADDRESS FOR THE POINTER TO THE NEXT GROUP USING VARIABLE 
C  IWK( LSVGRP + I ).
C
               IWK( LVALJR + I ) = J
               IWK( L )          = J + 1
  220      CONTINUE
         END IF
  230 CONTINUE
C
C  RESET THE STARTING ADDRESSES FOR THE LISTS OF GROUPS USING
C  EACH VARIABLE.
C
      DO 240 I    = N, 2, - 1
         L        = LSTAJC + I
         IWK( L ) = IWK( L - 1 )
  240 CONTINUE
      IWK( LSTAJC + 1 ) = 1
C
C  INITIALIZE WORKSPACE VALUES FOR SUBROUTINE HESPRD.
C
CDIR$ IVDEP
      DO 250 J             = 1, MAX( N, NG )
         IWK( LIUSED + J ) = 0
  250 CONTINUE
C
C  DEFINE FURTHER PARTITIONS OF THE WORKSPACE WHENEVER FINITE-
C  -DIFFERENCE GRADIENTS ARE USED.
C
      IF (  FDGRAD ) THEN
C
C  -- CONTINUE SETTING STARTING ADDRESSES FOR PARTITIONS OF IWK. --
C
C  THE RANGE TRANSFORMATION FOR EACH NONLINEAR ELEMENT IS OF A GIVEN
C  TYPE. SUPPOSE THERE ARE NTYPE NON-TRIVIAL TYPES. IWK( LSTYPE + I )
C  GIVES THE TYPE OF NONLINEAR ELEMENT I FOR I = 1, ...., NEL.
C
         LSTYPE = LSEND
C
C  THE RANGE TRANSFORMATION FROM ELEMENTAL TO INTERNAL VARIABLES IS
C  DEFINED BY A MATRIX W. FOR EACH NON-TRIVIAL TRANSFORMATION, THE
C  MATRIX W IS RECORDED. THE INFORMATION FOR THE I-TH TYPE STARTS IN
C  LOCATION LSWTRA + IWK( LSSWTR + I ), I = 1, ...., NTYPE.
C
         LSSWTR = LSTYPE + NEL
C
C  FOR EACH TYPE OF NONLINEAR ELEMENT USING A NONTRIVIAL RANGE 
C  TRANSFORMATION, INTEGER INFORMATION ISD ALSO RECORDED. 
C  THE INFORMATION FOR THE I-TH TYPE STARTS IN LOCATION
C  LSIWTR + IWK( LSSIWT + I ), I = 1, ...., NTYPE.
C
         LSSIWT = LSSWTR + NEL
C
C  THE FOLLOWING PIECES OF INTEGER INFORMATION ARE RECORDED ABOUT
C  THE I-TH TYPE OF NONLINEAR ELEMENT:
C
C    IWK( LSIWTR + IWK( LSSIWT + I ) + 1 ):
C            THE NUMBER OF INTERNAL VARIABLES, NINVR.
C    IWK( LSIWTR + IWK( LSSIWT + I ) + 2 ):
C            THE NUMBER OF ELEMENTAL VARIABLES, NELVR.
C    IWK( LSIWTR + IWK( LSSIWT + I ) + 2 + I ),
C         I = 1, ..., NELVR +  NINVR:
C            PIVOT SEQUENCES FOR THE LU FACTORS OF W.
C
C  AFTER THE FACTORIZATION AND COMPRESSION, ONLY NINVR LINEARLY
C  INDEPENDENT COLUMNS OF W ARE STORED.
C
         LSIWTR = LSSIWT + NEL
         IF ( LIWK .LT. LSIWTR ) THEN
            WRITE( IOUT, 2030 ) LSIWTR - LIWK
            INFORM = 4
            RETURN
         END IF
C 
C  -- CONTINUE SETTING STARTING ADDRESSES FOR PARTITIONS OF WK. --
C
C  THE FOLLOWING PIECES OF INTEGER INFORMATION ARE RECORDED ABOUT
C  THE I-TH TYPE OF NONLINEAR ELEMENT:
C
C    IWK( LSWTRA + IWK( LSSWTR + I ) + I ),
C         I = 1, ..., NELVR * NINVR:  THE MATRIX W STORED BY COLUMNS.
C
C  AFTER THE FACTORIZATION AND COMPRESSION, ONLY NINVR LINEARLY
C  INDEPENDENT COLUMNS OF W ARE STORED.
C
         LSWTRA = LWKSTR
C
C  ---------------------------------------------------------------------
C  CONSIDER ONLY ELEMENTS WHICH USE INTERNAL VARIABLES.
C  ---------------------------------------------------------------------
C
         NTYPE  = 0 
         LWFREE = 1
         LIWFRE = 1
C
C  LOOP OVER ALL NONLINEAR ELEMENTS.
C
         DO 350 IEL = 1, NEL
            IF ( INTREP( IEL ) ) THEN
C
C  CALCULATE THE RANGE TRANSFORMATION MATRIX W.
C 
               IS    = ISTAEV( IEL )
               NINVR = INTVAR( IEL + 1 ) - INTVAR( IEL )
               NELVR = ISTAEV( IEL + 1 ) - IS
               LW1   = LWFREE + NINVR * NELVR - 1
               L     = LSWTRA + LW1
CDIR$ IVDEP
               DO 280 I = 1, NELVR
                  WK( L + I ) = ZERO
  280          CONTINUE   
               K  = LSWTRA + LWFREE
               IS = IS - 1
               DO 320 I = 1, NELVR
                  WK( L + I ) = ONE
C                 write(6,*) ' L + 1, K ', L + 1, K 
                  CALL RANGES( IEL, .FALSE., WK( L + 1 ), 
     *                         WK( K ), NELVR, NINVR )
C                 WRITE(6,*) ' col ', I, ( WK( K + II - 1 ),
C    *                                     II = 1, NINVR )
                  WK( L + I ) = ZERO
                  K           = K + NINVR
C
C  CHECK TO SEE IF ANY OF THE COLUMNS BELONG TO DUPLICATED VARIABLES.
C
                  II = IELVAR( IS + I )
                  DO 290 J = 1, I - 1
                     IF ( IELVAR( IS + J ) .EQ. II ) GO TO 300
  290             CONTINUE
                  GO TO 320   
C
C  AMALGAMATE COLUMNS FROM DUPLICATE VARIABLES.
C
  300             CONTINUE
C                 write(6,*) 'AMALGAMATING VARIABLES ', I, J
                  KK = LSWTRA + LWFREE + ( J - 1 ) * NINVR - 1
                  LL = K - NINVR - 1
                  DO 310 JJ = 1, NINVR
C                       write(6,*) 'before wK(KK), wk(LL)',
C    *                              WK( KK + JJ ),WK( LL + JJ )
                     WK( KK + JJ ) = WK( KK + JJ ) + WK( LL + JJ )
                     WK( LL + JJ ) = ZERO
C                       write(6,*) 'after wK(KK), wk(LL)',
C    *                              WK( KK + JJ ),WK( LL + JJ )
  310             CONTINUE   
  320          CONTINUE
C
C  COMPARE THIS TRANSFORMATION MATRIX WITH PREVIOUS ONES.
C
               DO 340 I = 1, NTYPE
                  IF ( IWK( LSIWTR + IWK( LSSIWT + I ) ) .NE. NINVR 
     *                 .OR. IWK( LSIWTR + IWK( LSSIWT + I ) + 1 )
     *                 .NE. NELVR ) GO TO 340
                  DO 330 J = 0, NINVR * NELVR - 1
                     IF ( WK( LSWTRA + LWFREE + J ) .NE. WK( 
     *                    LSWTRA + IWK( LSSWTR + I ) + J ) ) GO TO 340
  330             CONTINUE
C
C  THE TRANSFORMATION IS AN EXISTING ONE. RECORD WHICH ONE.
C
                  IWK( LSTYPE + IEL ) = I
                  GO TO 350   
  340          CONTINUE
C
C  ENSURE THAT THERE IS SUFFICIENT ROOM.
C
               IF ( LIWFRE + 2 + NINVR + NELVR .GT. LIWK ) THEN
                  INFORM = 4
                  WRITE( IOUT, 2030 ) LIWFRE + 2 + NINVR + NELVR - LIWK 
                  RETURN
               END IF
               IF ( LWFREE + NINVR * NELVR .GT. LWK ) THEN
                  INFORM = 5
                  WRITE( IOUT, 2040 ) LWFREE + NINVR * NELVR - LWK 
                  RETURN
               END IF
C
C  THE TRANSFORMATION DEFINES A NEW TYPE. RECORD ITS DETAILS.
C
               NTYPE                      = NTYPE + 1
               IWK( LSTYPE + IEL        ) = NTYPE
               IWK( LSSIWT + NTYPE      ) = LIWFRE
               IWK( LSIWTR + LIWFRE     ) = NINVR
               IWK( LSIWTR + LIWFRE + 1 ) = NELVR
               IWK( LSSWTR + NTYPE      ) = LWFREE
               LIWFRE = LIWFRE + 2 + NINVR + NELVR
               LWFREE = LWFREE     + NINVR * NELVR 
            ELSE
              IWK( LSTYPE + IEL ) = 0
            END IF
  350    CONTINUE
C
C  FOR EACH TYPE OF ELEMENT WITH INTERNAL VARIABLES:
C
         DO 360 I  = 1, NTYPE
            LIWFRE = IWK( LSSIWT + I      )
            LWFREE = IWK( LSSWTR + I      )
            NINVR = IWK( LSIWTR + LIWFRE )
C
C  FACTORIZE W. USE GAUSSIAN ELIMINATION WITH COMPLETE PIVOTING.
C  DETERMINE THE "MOST INDEPENDENT" SET OF COLUMNS OF W.
C
CS          CALL SGELIM( NINVR, IWK( LSIWTR + LIWFRE + 1 ),
CD          CALL DGELIM( NINVR, IWK( LSIWTR + LIWFRE + 1 ),
     *                   IWK( LSIWTR + LIWFRE + 2 ),
     *                   IWK( LSIWTR + LIWFRE + NINVR + 2 ),
     *                   WK ( LSWTRA + LWFREE ) )
  360    CONTINUE   
C
C  COMPRESS THE DATA STRUCTURES TO REMOVE REDUNDANT INFORMATION.
C
         IF ( NTYPE .LT. NEL ) THEN
            K = LSSWTR + NTYPE
C 
C  COMPRESS INTEGER DATA.
C 
            DO 370 I = 1, NTYPE
               IWK( K + I ) = IWK( LSSIWT + I )
  370       CONTINUE
            LSSIWT = K
         END IF
         K      = LSSIWT + NTYPE
         LNIWTR = 0
         LNWTRA = 0
         DO 400 I = 1, NTYPE
            LIWFRO = IWK( LSSIWT + I      ) - 1
            NINVR = IWK( LSIWTR + LIWFRO + 1 )
            DO 380 J = 1, 2 * NINVR + 2
               IWK( K + LNIWTR + J ) = IWK( LSIWTR + LIWFRO + J )
  380       CONTINUE
            IWK( LSSIWT + I ) = LNIWTR + 1
            LNIWTR = LNIWTR + 2 + 2 * NINVR
C
C  COMPRESS REAL DATA.
C
            LWFREO   = IWK( LSSWTR + I ) - 1
            DO 390 J = 1, NINVR * NINVR
               WK( LSWTRA + LNWTRA + J ) = WK( LSWTRA + LWFREO + J )
  390       CONTINUE
            IWK( LSSWTR + I ) = LNWTRA + 1
            LNWTRA = LNWTRA + NINVR * NINVR 
  400    CONTINUE   
C
C  RECORD THE LENGTHS OF THE PARTITIONS OF THE WORKSPACE USED.
C
         LSIWTR = K
         LWKSTR = LSWTRA + LNWTRA
C
C  ---------------------------------------------------------------------
C  THE LIST OF VARIABLES IS ALLOCATED TO NSETS DISJOINTS SETS.
C  VARIABLE I OCCURS IN SET IWK( LSISET + I ).
C  ---------------------------------------------------------------------
C
         LSISET = LSIWTR + LNIWTR
         LSSVSE = LSISET + N
         LSEND  = LSSVSE + N + 1
         IF ( LIWK .LT. LSEND + N ) THEN
            WRITE( IOUT, 2030 ) LSEND + N - LIWK
            INFORM = 4
            RETURN
         END IF
C
C  ASSIGN INITIAL SET NUMBERS TO EACH VARIABLE.
C
         NSETS = 0
CDIR$    IVDEP
         DO 410 I = 1, N
            IWK( LSISET + I ) = N
  410    CONTINUE   
C
C  USE THE CURTIS-POWELL-REID ALGORITHM TO DETERMINE WHICH SET EACH
C  VARIABLE BELONGS TO. LOOP OVER THE VARIABLES.
C
         DO 500 I = 1, N
C
C  LOOP OVER THE ELEMENTS WHICH USE VARIABLE I.
C  THE ELEMENTS ARE OBTAINED FROM A LINKED-LIST.
C
            VRUSED = .FALSE.
            IPT    = IWK( LSPTRS + I )
            IF ( IPT .GE. 0 ) THEN
               IELL = IWK( LSELTS + I ) 
  420          CONTINUE
               IEL   = IELING( IELL )
               ITYPE = IWK( LSTYPE + IEL )
C              WRITE( 6, * ) ' ELEMENT ', IEL
C
C  CHECK THAT THE VARIABLE BELONGS TO THE "INDEPENDENCE" SET OF
C  ELEMENTS WITH INTERNAL VARIABLES.
C
               IF ( ITYPE .GT. 0 ) THEN
                  LIWFRE = IWK( LSSIWT + ITYPE  )
                  NINVR  = IWK( LSIWTR + LIWFRE )
                  DO 430 J = 1, NINVR
                     K = J - 1
                     L = IWK( LSIWTR + LIWFRE + NINVR + 1 + J ) - 1
                     IF ( I .EQ. IELVAR( ISTAEV( IEL ) + L ) )
     *                  GO TO 440
  430             CONTINUE        
                  GO TO 450
  440             CONTINUE        
               END IF 
               VRUSED = .TRUE.
  450          CONTINUE
C
C  LOOP OVER THE COMPLETE LIST OF VARIABLES USED BY ELEMENT IEL.
C
CDIR$ IVDEP
               DO 460 J = ISTAEV( IEL ), ISTAEV( IEL + 1 ) - 1
C 
C  IF VARIABLE IV IS USED, FLAG THE SET THAT CONTAINS IT.
C
                  IWK( LSEND + IWK( LSISET + IELVAR( J ) ) ) = 1
  460          CONTINUE
C
C  CHECK THE LINK-LIST TO SEE IF FURTHER ELEMENTS USE THE VARIABLE.
C
               IF ( IPT .GT. 0 ) THEN
                  IELL = IWK( LSELTS + IPT )
                  IPT  = IWK( LSPTRS + IPT )
                  GO TO 420
               END IF
            END IF 
C
C  SEE IF THE VARIABLE MAY BE PLACED IN THE FIRST NSETS SETS.
C
            IF ( VRUSED ) THEN
               DO 470 J = 1, NSETS
                  IF ( IWK( LSEND + J ) .EQ. 0 ) GO TO 480
                  IWK( LSEND + J ) = 0
  470          CONTINUE
C
C  THE VARIABLE NEEDS A NEW SET.
C
               NSETS = NSETS + 1
               J     = NSETS
C
C  THE VARIABLE WILL BE PLACED IN SET J.
C
  480          CONTINUE
               IWK( LSISET + I ) = J        
C
C  RESET THE FLAGS TO ZERO.
C
CDIR$ IVDEP
               DO 490 K = J, NSETS
                  IWK( LSEND + K ) = 0
  490          CONTINUE
            ELSE
C
C  THE VARIABLE IS NOT TO BE USED.
C
               IWK( LSISET + I ) =  N
            END IF
  500    CONTINUE   
C
C  CHECK THAT THERE IS AT LEAST ONE SET.
C
         IF ( NSETS .NE. 0 ) THEN
C
C  PRINT OUTPUT.
C
CDIR$ IVDEP
            DO 510 I = 1, N
               IWK( LSISET + I ) = MIN( IWK( LSISET + I ), NSETS + 1 )
C              WRITE( 6, * ) ' VARIABLE ', I, ' SET ', 
C    *                       IWK( LSISET + I )
  510       CONTINUE
C
C  ---------------------------------------------------------------------
C  OBTAIN A LIST, IWK(LSEND), OF THE VARIABLES CORRESPONDING TO EACH SET.
C  ---------------------------------------------------------------------
C
C  CLEAR IWK( LSSVSE ).
C
CDIR$ IVDEP
            DO 520 J = 2, NSETS + 2
              IWK( LSSVSE + J ) = 0
  520       CONTINUE
C
C  COUNT THE NUMBER OF ELEMENTS IN EACH SET AND STORE IN IWK( LSSVSE ).
C  NEGATE THE SET NUMBERS IN IWK( LSISET ), SO THAT THEY ARE FLAGGED
C  AS IWK( LSISET ) IS GRADUALLY OVERWRITTEN BY VARIABLE INDICES.
C
            DO 530 K = 1, N
               J     = IWK( LSISET + K )
               IWK( LSISET + K )     = - J
               IWK( LSSVSE + J + 1 ) = IWK( LSSVSE + J + 1 ) + 1
  530       CONTINUE
C
C  COMPUTE THE STARTING ADDRESSES FOR EACH SET WITHIN IWK( LSISET ).
C
            IWK( LSSVSE + 1 ) = 1
            DO 540 J          = 2, NSETS + 2
               IWK( LSSVSE + J ) = IWK( LSSVSE + J ) +
     *                             IWK( LSSVSE + J - 1 )
  540       CONTINUE
C
C  STORE IN IWK( LSEND ) THE VARIABLE WHOSE SET NUMBER
C  IS THE IWK( LSSVSE + J )-TH ENTRY OF IWK( LSEND ).
C
            ISOFAR   = 0
            DO 570 J = 1, NSETS + 1
               ISTARJ = IWK( LSSVSE + J )
               DO 550 IVARP1 = ISOFAR + 1, N
                  IF ( ISTARJ .LT. IVARP1 ) GO TO 560
  550          CONTINUE
               IVARP1 = N + 1
  560          CONTINUE
               ISOFAR = IVARP1 - 1
               IWK( LSEND + J ) = ISOFAR
  570       CONTINUE
C
C  REORDER THE ELEMENTS INTO SET ORDER.
C  FILL IN EACH SET FROM THE FRONT. AS A NEW ENTRY IS PLACED
C  IN SET K INCREASE THE POINTER IWK( LSSVSE + K ) BY ONE AND FIND
C  THE NEW VARIABLE, IWK( LSEND + K ), THAT CORRESPONDS TO THE SET NOW
C  POINTED TO BY IWK( LSSVSE + K ).
C
            DO 660 J = 1, NSETS + 1
C
C  DETERMINE THE NEXT UNPLACED ENTRY, IWK( LSSVSE ), IN IWK( LSISET ).
C
  610          CONTINUE
               ISTRT = IWK( LSSVSE + J )
C
C  SEE IF ALL THE ELEMENTS IN SET J HAVE BEEN ASSIGNED.
C  
               IF ( ISTRT .EQ. IWK( LSSVSE + J + 1 ) ) GO TO 660
               IF ( IWK( LSISET + ISTRT ) .GT. 0 ) GO TO 660
C
C  EXTRACT THE VARIABLE AND SET NUMBERS OF THE STARTING ELEMENT.
C
               IVAR = IWK( LSEND + J )
               JSET = - IWK( LSISET + ISTRT )
C
C  MOVE ELEMENTS IN A CYCLE, ENDING BACK AT SET J.
C
               DO 640 K = ISTRT, N
C
C  FIND THE FIRST EMPTY LOCATION IN SET JSET IN IWK( LSEND )
C
                  INEXT = IWK( LSSVSE + JSET )
C
C  EXTRACT THE VARIABLE INDEX OF THE NEXT ELEMENT.
C
                  NEWVAR = IWK( LSEND + JSET )
C
C  UPDATE IWK( LSSVSE + JSET ), FIND THE NEW VARIABLE INDEX AND STORE
C  IT IN IWK( LSEND + JSET ).
C
                  ISTARJ         = INEXT + 1
                  IWK( LSSVSE + JSET ) = ISTARJ
                  DO 620 IVARP1  = NEWVAR + 1, N
                     IF ( ISTARJ .LT. IVARP1 ) GO TO 630
  620             CONTINUE
                  IVARP1 = N + 1
  630             CONTINUE
                  IWK( LSEND + JSET ) = IVARP1 - 1
C
C  IF THE ENTRY BELONGS IN THE J-TH SET, THE CYCLE IS COMPLETE.
C
                  IF ( JSET .EQ. J ) GO TO 650
C
C  EXTRACT THE NUMBER OF THE SET OF THE NEXT ELEMENT.
C
                  NEWSET = - IWK( LSISET + INEXT )
C
C  STORE THE VARIABLE INDEX OF THE CURRENT ELEMENT.
C
                  IWK( LSISET + INEXT ) = IVAR
C
C  MAKE THE NEXT ELEMENT INTO THE CURRENT ONE.
C
                  IVAR = NEWVAR
                  JSET = NEWSET
  640          CONTINUE
C
C  THE CYCLE IS COMPLETE.
C
  650          CONTINUE
C
C  STORE THE VARIABLE INDEX OF THE STARTING ELEMENT.
C
               IWK( LSISET + ISTRT ) = IVAR
               GO TO 610
  660       CONTINUE
C
C  REVISE IWK( LSSVSE ) TO POINT TO THE START OF EACH SET.
C
            DO 670 J             = NSETS + 1, 2, - 1
               IWK( LSSVSE + J ) = IWK( LSSVSE + J - 1 )
  670       CONTINUE
            IWK( LSSVSE + 1 ) = 1
         END IF
C        DO 671 I = 1, NSETS
C           WRITE( 6, * ) ' SET ', I, ' VARIABLES ', 
C    *       ( IWK( LSISET + J ), J = IWK( LSSVSE + I ), 
C    *         IWK( LSSVSE + I + 1 ) - 1 )
C 671    CONTINUE   
      ELSE
C
C  EXACT GRADIENTS ARE USED. NO FURTHER PARTITIONING OF THE WORKSPACE
C  IS NEEDED.
C
         LSTYPE = LSEND
         LSSWTR = LSTYPE
         LSSIWT = LSSWTR
         LSIWTR = LSSIWT
         LSISET = LSIWTR
         LSSVSE = LSISET
         LSWTRA = LWKSTR
      END IF    
C
C  SET THE LENGTH OF THE REMAINING PARTITIONS OF THE WORKSPACE FOR
C  ARRAY BOUND CHECKING IN CALLS TO OTHER SUBPROGRAMS.
C
      LNTYPE = MAX( 1, LSSWTR - LSTYPE )
      LNSWTR = MAX( 1, LSSIWT - LSSWTR )
      LNSIWT = MAX( 1, LSIWTR - LSSIWT )
      LNIWTR = MAX( 1, LSISET - LSIWTR )
      LNISET = MAX( 1, LSSVSE - LSISET )
      LNSVSE = MAX( 1, LSEND  - LSSVSE )
      LNWTRA = MAX( 1, LWKSTR - LSWTRA )
C
C  RECORD THE LENGTHS OF THE REMAINING INTEGER AND REAL WORKSPACE.
C
      LIWK2  = LIWK - LSEND
      LWK2   = LWK  - LWKSTR
C
C  -- SET THE STARTING ADDRESSES FOR THE PARTITIONS WITHIN FUVALS. --
C
C  A FULL DESCRIPTION OF THE PARTITIONS OF FUVALS IS GIVEN IN THE 
C  THE INTRODUCTORY COMMENTS TO SUBROUTINE SBMIN.
C
      LFXI   = 0
      LGXI   = LFXI           + NEL
      LHXI   = INTVAR( NEL1 ) - 1
      LGGFX  = LGGFX          - 1
      LDX    = LGGFX          + N
      LGRJAC = LDX            + N
      LEND   = LGRJAC         + NVARGP
C
C  PRINT ALL OF THE STARTING ADDRESSES FOR THE WORKSPACE ARRAY 
C  PARTITIONS.
C
      IF ( IPRINT .GE. 3 ) THEN
         WRITE( IOUT, 2000 )
     *     LFXI, LGXI, LHXI, LGGFX, LDX, LGRJAC, LEND, LFUVAL
         WRITE( IOUT, 2010 ) LSPTRS, LSELTS, LINDEX,
     *                       LSWKSP, LIUSED, LFREEC, LNNONZ, LNONZ2,
     *                       LSYMMD, LSYMMH, LSLGRP, LSTAJC,
     *                       LSTAGV, LSVGRP, LGCOLJ, LVALJR, 
     *                       LSTYPE, LSSWTR, LSSIWT, LSIWTR, LSISET,
     *                       LSSVSE, LSEND , LIWK 
         WRITE( IOUT, 2020 ) LQGRAD, LBREAK, LP,     LXCP, LX0, LGX0, 
     *                       LDELTX, LBND  , LSWTRA, LWKSTR, LWK
      END IF
C
C  CHECK THAT THE ARRAY FUVALS HAS SUFFICIENT ROOM FOR THE CALCULATION.
C
      IF ( LEND .GT. LFUVAL ) THEN
         WRITE( IOUT, 2050 ) LEND - LFUVAL
         INFORM = 5
         RETURN
      END IF
C
C  SET THE LENGTH OF EACH PARTITION OF THE REAL WORKSPACE ARRAY 
C  FUVALS FOR ARRAY BOUND CHECKING IN CALLS TO OTHER SUBPROGRAMS.
C
      LNFXI  = MAX( 1, LGXI   - LFXI   )
      LNGXI  = MAX( 1, LHXI   - LGXI   )
      LNGUVL = MAX( 1, LHXI   - LFXI   )
      LNHXI  = MAX( 1, LGGFX  - LHXI   )
      LNHUVL = MAX( 1, LGGFX  - LFXI   )
      LNGGFX = MAX( 1, LDX    - LGGFX  ) 
      LNDX   = MAX( 1, LGRJAC - LDX    )
      LNGRJC = MAX( 1, LEND   - LGRJAC )
C ** Correction 1. 05/02/93: 2 lines added ** 
      MAXSIN = MAX( 1, MAXSIN )
      MAXSEL = MAX( 1, MAXSEL )
C ** Correction 1. 05/02/93: end of correction **
      INFORM = 0
C
C  NON-EXECUTABLE STATEMENTS.
C
 2000 FORMAT( /, ' Starting addresses for the partitions of FUVALS ',
     *        /, ' ----------------------------------------------- ',
     *       //, '   LFXI   LGXI   LHXI  LGGFX ',
     *           '   LDX LGRJAC   LEND   LFUVAL ', /, 7I7, I9 )
 2010 FORMAT( /, ' Starting addresses for partitions of IWK ',
     *        /, ' ---------------------------------------- ', //,
C ** Correction 2. 18/07/94: 1 line corrected ** 
     *  ' LSPTRS LSELTS LINDEX LSWKSP LIUSED LFREEC LNNONZ LNONZ2  ...',
C ** Correction 2. 18/07/94: end of correction **
     *   /, 8I7, //,
     *   ' ...... LSYMMD LSYMMH LSLGRP LSTAJC LSTAGV LSVGRP LGCOLJ ...',
     *   /, 7X, 7I7, //,
C ** Correction 3. 18/07/94: 1 line corrected ** 
     *  ' ...... LVALJR LSTYPE LSSWTR LSSIWT LSIWTR LSISET LSSVSE ... ',
C ** Correction 3. 18/07/94: end of correction **
     *   /, 7X, 7I7, //,
     *   ' ......  LSEND     LIWK ',
     *   /, 7X, I7, I9 )
 2020 FORMAT( /, ' Starting addresses for partitions of WK ',
     *        /, ' --------------------------------------- ', //,
     *   '   LQGRAD   LBREAK       LP     LXCP      LX0     LGX0 ... ',
     *   /, 6I9, //,
     *   '   ......   LDELTX     LBND   LSWTRA   LWKSTR      LWK ',
     *   /, 9X, 5I9 )
 2030 FORMAT( /, ' The size of array IWK must be increased',  
     *           ' by at least ', I8 )
 2040 FORMAT( /, ' The size of array WK must be increased',  
     *           ' by at least ', I8 )
 2050 FORMAT( /, ' The size of array FUVALS must be increased',  
     *           ' by at least ', I8 )
      RETURN
C
C  END OF SUBROUTINE INITW.
C
      END
