* SUBROUTINE PF1HS2                ALL SYSTEMS                99/12/01
C PORTABILITY : ALL SYSTEMS
C 99/12/01 TU : ORIGINAL VERSION
*
* PURPOSE :
* NUMERICAL COMPUTATION OF THE HESSIAN MATRIX OF THE MODEL FUNCTION
* USING ITS GRADIENTS - SPARSE VERSION USING DIRECT COLOURING METHOD.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  II  ML SIZE OF THE COMPACT FACTOR.
*  II  M  NUMBER OF NONZERO ELEMENTS OF THE SPARSE HESSIAN MATRIX.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  RA  XO(NF)  AUXILIARY VECTOR.
*  RO  HF(M)  HESSIAN MATRIX OF THE MODEL FUNCTION.
*  IU  IH(NF+1)  POINTER VECTOR OF SPARSE HESSIAN MATRIX.
*  IU  JH(M)  INDEX VECTOR OF THE HESSIAN MATRIX.
*  RI  GF(NF)  GRADIENT OF THE MODEL FUNCTION.
*  RA  GO(NF)  AUXILIARY VECTOR.
*  II  COL(NF)  VECTOR DISCERNING GROUPS OF THE HESSIAN COLUMN OF THE
*         SAME COLOUR.
*  IA  WN11(NF+1)  AUXILIARY VECTOR.
*  IA  WN12(NF+1)  AUXILIARY VECTOR.
*  RA  XS(NF)  AUXILIARY VECTOR USED FOR STEP SIZES.
*  RI  FF  VALUE OF THE MODEL FUNCTION.
*  RI  ETA1  PRECISION OF THE COMPUTED VALUES.
*  II  KBF  TYPE OF BOUNDS. KBF=0-BOUNDS ARE NOT USED. KBF=1-ONE SIDED
*         BOUNDS. KBF=1-TWO SIDED BOUNDS.
*  IU  ITERM  TERMINATION INDICATOR.
*  IU  ISYS  CONTROL PARAMETER.
*
* SUBPROGRAMS USED :
*  S   MXSTG1  WIDTHEN THE STRUCTURE.
*  S   MXSTL2  SHRINK THE STRUCTURE.
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVSET  INITIATION OF A VECTOR.
*
      SUBROUTINE PF1HS2(NF,ML,M,X,IX,XO,HF,IH,JH,GF,GO,COL,WN11,
     & WN12,XS,FF,ETA1,KBF,ITERM,ISYS)
      INTEGER NF,ML,M,IX(NF),IH(NF+1),JH(M),COL(NF),WN11(NF+1),
     & WN12(NF+1),KBF,ITERM,ISYS
      DOUBLE PRECISION X(NF),XO(NF),HF(M),GF(NF),GO(NF),XS(NF),
     & FF,ETA1
      DOUBLE PRECISION XTEMP,FTEMP,ETA
      INTEGER I,J,J1,K,K1,L,MX,MM,IVAR,JVAR
      SAVE MX,MM,IVAR,JVAR
      SAVE XTEMP,FTEMP,ETA
      IF(ITERM.NE.0) GO TO 12
      IF (ISYS.EQ.1) GO TO 3
C      GO TO (1,3) ISYS+1
    1 CONTINUE
      MM=IH(NF+1)-1
      IF(3*MM-NF+ML.GE.M) THEN
      ITERM=-45
      ISYS=0
      RETURN
      END IF
      ETA=SQRT(ETA1)
      FTEMP=FF
      CALL MXVCOP(NF,X,XO)
C
C     WIDTHEN THE STRUCTURE
C
      K=2*MM-NF
      DO 50 I=ML+MM,1,-1
        JH(K+I)=JH(MM+I)
   50 CONTINUE
      CALL MXSTG1(NF,MX,IH,JH,WN12,WN11)
      CALL MXVSET(K,0.0D 0,HF)
      IVAR=1
    2 CONTINUE
      IF(IVAR.GT.NF) GO TO 870
      DO 200 J=IVAR,NF
        IF(COL(J).GE.1) THEN
          GO TO 200
        ELSE
          JVAR=J
          GO TO 300
        END IF
 200  CONTINUE
 300  CONTINUE
      DO 400 J=IVAR,JVAR
        L=ABS(COL(J))
        IF(KBF.GT.0) THEN
          IF(IX(L).LE.-7) GO TO 400
        END IF
C     STEP SELECTION
C
        XS(L)=ETA*MAX(ABS(X(L)),1.0D 0)*SIGN(1.0D 0,X(L))
        XTEMP=X(L)
        X(L)=XTEMP+XS(L)
        XS(L)=X(L)-XTEMP
 400  CONTINUE
        ISYS=1
        RETURN
    3 CONTINUE
C
C     NUMERICAL DIFFERENTIATION
C
C
C     SET AUXILIARY VECTOR DISCERNING THE SINGLETONS IN A GROUP TO ZERO
C
      DO 450 J1=1,NF
        WN11(J1)=0
  450 CONTINUE
C
C     DISCERN SINGLETONS OF THE GROUP OF THE SAME COLOR.
C
      DO 600 J1=IVAR,JVAR
        L=ABS(COL(J1))
        DO 500 K=IH(L),IH(L+1)-1
          K1=ABS(JH(K))
          IF(WN11(K1).EQ.0) THEN
            WN11(K1)=J1
          ELSE
            WN11(K1)=-1
          END IF
 500  CONTINUE
 600  CONTINUE
C
C     NUMERICAL VALUES COMPUTATION
C
      DO 800 J1=IVAR,JVAR
        L=ABS(COL(J1))
        DO 700 K=IH(L),IH(L+1)-1
          K1=ABS(JH(K))
          IF(WN11(K1).GT.0) THEN
            HF(K)=(GF(K1)-GO(K1))/XS(L)
          END IF
 700    CONTINUE
 800  CONTINUE
C
C     SET THE ORIGINAL VALUE OF X FOR THE COMPONENTS OF THE ACTUAL COLOR.
C
      DO 850 J=IVAR,JVAR
        L=ABS(COL(J))
        X(L)=XO(L)
 850  CONTINUE
      IVAR=JVAR+1
      GO TO 2
 870  CONTINUE
C
C     MOVE THE ELEMENTS OF THE HESSIAN APPROXIMATION INTO THE UPPER
C     TRIANGULAR PART
C
      DO 900 I=1,NF
        WN11(I)=WN12(I)+1
 900  CONTINUE
      DO 1100 I=1,NF
        IVAR=IH(I)
        JVAR=WN12(I)-1
        DO 1000 J=IVAR,JVAR
          K=ABS(JH(J))
          L=WN11(K)
          IF(HF(L).EQ.0) THEN
            HF(L)=HF(J)
          ELSE IF(HF(L).NE.0.AND.HF(J).NE.0) THEN
            HF(L)=0.5D 0*(HF(J)+HF(L))
          END IF
          WN11(K)=WN11(K)+1
 1000   CONTINUE
 1100 CONTINUE
    5 FF=FTEMP
C
C     SHRINK THE STRUCTURE
C
      CALL MXSTL2(NF,MX,HF,IH,JH,WN12)
      K=2*MM-NF
      DO 1200 I=1,ML+MM
        JH(MM+I)=JH(K+I)
 1200 CONTINUE
C
C     RETRIEVE VALUES
C
      CALL MXVCOP(NF,XO,X)
   12 CONTINUE
      ISYS=0
      RETURN
      END
* SUBROUTINE PFSED7             ALL SYSTEMS                   98/12/01
C PORTABILITY : ALL SYSTEMS
C 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* COMPUTATION OF THE SPARSE SCHUR COMPLEMENT STRUCTURE FOR THE
* KARUSH-KUHN-TUCKER MATRIX WITH THE UNIT HESSIAN BLOCK.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  II  NC  NUMBER OF CONSTRAINT FUNCTIONS.
*  IO  ND  ACTUAL NUMBER OF DENSE ROWS.
*  IO  MD  MAXIMUM NUMBER OF DENSE ROWS.
*  IO  MDE  MINIMUM NUMBER OF NONZERO ELEMENTS IN DENSE ROWS.
*  IO  IB(NF+1)  POINTERS OF DIAGONAL ELEMENTS OF THE SPARSE SCHUR
*         COMPLEMENT.
*  IO  JB(M)  COLUMN INDICES OF NONZERO ELEMENTS OF THE SPARSE SCHUR
*         COMPLEMENT.
*  II  ICG(NC+1)  POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD CG.
*  II  JCG(MC)  COLUMN INDICES OF ELEMENTS IN THE FIELD CG.
*  IO  KCG(NF)  INDICATION OF DENSE ROWS.
*  IO  ITERM  TERMINATION INDICATOR.
*
      SUBROUTINE PFSED7(NF,NC,ND,MD,MDE,IB,JB,ICG,JCG,KCG,ITERM)
      INTEGER NF,NC,ND,MD,MDE,IB(*),JB(*),ICG(*),JCG(*),KCG(*)
      INTEGER K,KK,K1,K2,L,LL,L2,MM,KC,KCP,ITERM
      CALL MXVINS(NF,0,KCG)
      K2=0
      DO 2 KC=1,NC
      K1=K2+1
      K2=ICG(KC+1)-1
      DO 1 KK=K1,K2
      K=ABS(JCG(KK))
      KCG(K)=KCG(K)+1
    1 CONTINUE
    2 CONTINUE
      ND=0
      DO 3 K=1,NF
      IF (KCG(K).GE.MDE) THEN
      ND=ND+1
      KCG(K)=ND
      ELSE
      KCG(K)=0
      END IF
    3 CONTINUE
      IF (ND.GT.MD) THEN
      ITERM=-46
      RETURN
      END IF
      MM=0
      IB(1)=1
      DO 6 KC=1,NC
      MM=MM+1
      JB(MM)=KC
      K1=ICG(KC)
      K2=ICG(KC+1)-1
      DO 5 KCP=KC+1,NC
      KK=K1
      LL=ICG(KCP)
      L2=ICG(KCP+1)-1
      K=ABS(JCG(KK))
      L=ABS(JCG(LL))
    4 IF (KCG(K).NE.0.OR.K.LT.L) THEN
      IF (KK.GE.K2) GO TO 5
      KK=KK+1
      K=ABS(JCG(KK))
      GO TO 4
      ELSE IF (KCG(L).NE.0.OR.L.LT.K) THEN
      IF (LL.GE.L2) GO TO 5
      LL=LL+1
      L=ABS(JCG(LL))
      GO TO 4
      ELSE
      MM=MM+1
      JB(MM)=KCP
      END IF
    5 CONTINUE
      IB(KC+1)=MM+1
    6 CONTINUE
      RETURN
      END
* SUBROUTINE PFSED9             ALL SYSTEMS                   98/12/01
C PORTABILITY : ALL SYSTEMS
C 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* COMPUTATION OF THE SPARSE SCHUR COMPLEMENT FOR THE KARUSH-KUHN-TUCKER
* MATRIX WITH THE DIAGONAL HESSIAN BLOCK.
*
* PARAMETERS :
*  II  NC  NUMBER OF CONSTRAINT FUNCTIONS.
*  II  ND  ACTUAL NUMBER OF DENSE ROWS.
*  RO  B(M)  ELEMENTS OF THE SPARSE MATRIX B.
*  II  IB(NC+1)  POINTERS OF THE DIAGONAL ELEMENTS OF B.
*  II  JB(M)  INDICES OF THE NONZERO ELEMENTS OF B.
*  II  CG(MC)  ELEMENTS OF THE CONSTRAINT JACOBIAN MATRIX.
*  II  ICG(NC+1)  POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD CG.
*  II  JCG(MC)  COLUMN INDICES OF ELEMENTS IN THE FIELD CG.
*  II  KCG(NF)  INDICATION OF DENSE ROWS.
*  RO  CGD(NC*ND)  ELEMENTS OF DENSE ROWS.
*  RI  DD(ND)  DENSE PART OF DIAGONAL WEIGHTING MATRIX.
*  RI  D(NF)  DIAGONAL WEIGHTING MATRIX.
*
      SUBROUTINE PFSED9(NC,ND,B,IB,JB,CG,ICG,JCG,KCG,CGD,DD,D)
      INTEGER NC,ND,IB(*),JB(*),ICG(*),JCG(*),KCG(*)
      DOUBLE PRECISION B(*),CG(*),CGD(*),DD(*),D(*)
      INTEGER I1,I2,II,K,KK,K1,K2,L,LL,L2,KC
      IF (ND.GT.0) CALL MXVSET(NC*ND,0.0D 0,CGD)
      DO 4 KC=1,NC
      I1=IB(KC)
      I2=IB(KC+1)-1
      K1=ICG(KC)
      K2=ICG(KC+1)-1
      DO 2 II=I1,I2
      B(II)=0.0D 0
      L=JB(II)
      KK=K1
      LL=ICG(L)
      L2=ICG(L+1)-1
      K=JCG(KK)
      L=JCG(LL)
    1 IF (KCG(K).NE.0.OR.K.LE.0.OR.K.LT.L) THEN
      IF (KK.GE.K2) GO TO 2
      KK=KK+1
      K=JCG(KK)
      GO TO 1
      ELSE IF (KCG(L).NE.0.OR.L.LE.0.OR.L.LT.K) THEN
      IF (LL.GE.L2) GO TO 2
      LL=LL+1
      L=JCG(LL)
      GO TO 1
      ELSE
      B(II)=B(II)+D(K)*CG(KK)*CG(LL)
      IF(KK.GE.K2.OR.LL.GE.L2) GO TO 2
      KK=KK+1
      LL=LL+1
      K=JCG(KK)
      L=JCG(LL)
      GO TO 1
      END IF
    2 CONTINUE
      IF (ND.GT.0) THEN
      DO 3 KK=K1,K2
      K=JCG(KK)
      LL=KCG(K)
      IF (LL.GT.0) THEN
      IF (K.GT.0) CGD((LL-1)*NC+KC)=CG(KK)
      DD(LL)=1.0D 0/D(K)
      END IF
    3 CONTINUE
      END IF
    4 CONTINUE
      RETURN
      END
* SUBROUTINE PFSET2             ALL SYSTEMS                   97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* COMPUTATION OF THE NUMBER OF NONZERO ELEMENTS OF THE SPARSE
* HESSIAN MATRIX STORED IN THE BLOCKED FORM.
*
* PARAMETERS :
*  II  NA  NUMBER OF APPROXIMATED FUNCTIONS.
*  II  MB  NUMBER OF NONZERO ELEMENTS OF THE PARTITIONED HESSIAN MATRIX
*  II  MC  MAXIMUM NUMBER OF ELEMENTS OF THE PARTIAL HESSIAN MATRIX.
*  RI  IAG(NA+1)  POSITION OF THE FIRST ROWS ELEMENTS IN THE SPARSE
*         JACOBIAN MATRIX.
*
      SUBROUTINE PFSET2(NA,MB,MC,IAG)
      INTEGER NA,MB,MC,IAG(*)
      INTEGER K,L,KA
      MB=0
      MC=0
      DO 1 KA=1,NA
      K=IAG(KA)
      L=IAG(KA+1)-K
      MB=MB+L*(L+1)/2
      MC=MAX(MC,L*(L+1)/2)
    1 CONTINUE
      RETURN
      END
* SUBROUTINE PFSET3             ALL SYSTEMS                   97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* COMPUTATION OF THE SPARSE STRUCTURE OF THE HESSIAN MATRIX FROM THE
* SPARSE STRUCTURE OF THE JACOBIAN MATRIX.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIADLES.
*  II  NA  NUMBER OF APPROXIMATED FUNCTIONS.
*  IO  IH(NF+1)  POINTERS OF THE DIAGONAL ELEMENTS OF H.
*  IO  JH(M)  INDICES OF THE NONZERO ELEMENTS OF H.
*  II  IAG(NA+1)  POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD AG.
*  II  JAG(MA)  COLUMN INDICES OF ELEMENTS IN THE FIELD AG.
*
      SUBROUTINE PFSET3(NF,NA,IH,JH,IAG,JAG)
      INTEGER NF,NA,IH(*),JH(*),IAG(*),JAG(*)
      INTEGER I,J,JF,JA,K,LF,LA,KA,M
      M=IH(NF+1)-1
      DO 7 KA=1,NA
      LA=IAG(KA+1)-1
      DO 6 K=IAG(KA),LA
      I=JAG(K)
      JF=IH(I)
      LF=IH(I+1)-1
      DO 5 JA=K,LA
      J=JAG(JA)
    2 IF (JH(JF).LT.J.AND.JF.LE.LF) THEN
      JF=JF+1
      IF (JF.LE.LF) GO TO 2
      END IF
      IF (JH(JF).GT.J .OR.JF.GT.LF) THEN
      DO 3 J=I+1,NF+1
      IH(J)=IH(J)+1
    3 CONTINUE
      DO 4 J=M,JF,-1
      JH(J+1)=JH(J)
    4 CONTINUE
      JH(JF)=JAG(JA)
      JF=JF+1
      LF=LF+1
      M=M+1
      END IF
    5 CONTINUE
    6 CONTINUE
    7 CONTINUE
      RETURN
      END
* SUBROUTINE PNINT3                ALL SYSTEMS                91/12/01
C PORTABILITY : ALL SYSTEMS
C 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* EXTRAPOLATION OR INTERPOLATION FOR LINE SEARCH WITHOUT DIRECTIONAL
* DERIVATIVES.
*
* PARAMETERS :
*  RI  RO  INITIAL VALUE OF THE STEPSIZE PARAMETER.
*  RI  RL  LOWER VALUE OF THE STEPSIZE PARAMETER.
*  RI  RU  UPPER VALUE OF THE STEPSIZE PARAMETER.
*  RI  RI  INNER VALUE OF THE STEPSIZE PARAMETER.
*  RI  FO  VALUE OF THE OBJECTIVE FUNCTION FOR R=RO.
*  RI  FL  VALUE OF THE OBJECTIVE FUNCTION FOR R=RL.
*  RI  FU  VALUE OF THE OBJECTIVE FUNCTION FOR R=RU.
*  RI  FI  VALUE OF THE OBJECTIVE FUNCTION FOR R=RI.
*  RO  PO  INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RO  R  VALUE OF THE STEPSIZE PARAMETER OBTAINED.
*  II  MODE  MODE OF LINE SEARCH.
*  II  MTYP  METHOD SELECTION. MTYP=1-BISECTION. MTYP=2-TWO POINT
*         QUADRATIC INTERPOLATION. MTYP=2-THREE POINT QUADRATIC
*         INTERPOLATION.
*  IO  MERR  ERROR INDICATOR. MERR=0 FOR NORMAL RETURN.
*
* METHOD :
* EXTRAPOLATION OR INTERPOLATION WITH STANDARD MODEL FUNCTIONS.
*
      SUBROUTINE PNINT3(RO,RL,RU,RI,FO,FL,FU,FI,PO,R,MODE,MTYP,MERR)
      DOUBLE PRECISION RO,RL,RU,RI,FO,FL,FU,FI,PO,R
      INTEGER MODE,MTYP,MERR,NTYP
      DOUBLE PRECISION AL,AU,AI,DEN,DIS
      LOGICAL L1,L2
      DOUBLE PRECISION ZERO,HALF,ONE,TWO,THREE,FOUR,C1L,C1U,C2L,C2U,C3L
      PARAMETER(ZERO=0.0D 0,HALF=0.5D 0,ONE=1.0D 0,TWO=2.0D 0,
     &  THREE=3.0D 0,FOUR=4.0D 0,C1L=1.1D 0,C1U=1.0D 3,
     &  C2L=1.0D-2,C2U=0.9D 0,C3L=1.0D-1)
      MERR = 0
      IF (MODE .LE. 0)  RETURN
      IF (PO .GE. ZERO)  THEN
      MERR = 2
      RETURN
      ELSE  IF (RU .LE. RL)  THEN
      MERR = 3
      RETURN
      END IF
      L1 = RL .LE. RO
      L2 = RI .LE. RL
      DO 1  NTYP = MTYP, 1, -1
      IF (NTYP .EQ. 1)  THEN
C
C     BISECTION
C
      IF (MODE .EQ. 1)  THEN
      R = TWO * RU
      RETURN
      ELSE IF (RI-RL.LE.RU-RI) THEN
      R=HALF*(RI+RU)
      RETURN
      ELSE
      R=HALF*(RL+RI)
      RETURN
      END IF
      ELSE IF (NTYP.EQ.MTYP.AND.L1) THEN
      IF (.NOT.L2) AI=(FI-FO)/(RI*PO)
      AU=(FU-FO)/(RU*PO)
      END IF
      IF (L1.AND.(NTYP.EQ.2.OR.L2)) THEN
C
C     TWO POINT QUADRATIC EXTRAPOLATION OR INTERPOLATION
C
      IF (AU.GE.ONE) GO TO 1
      R=HALF*RU/(ONE-AU)
      ELSE IF (.NOT.L1.OR..NOT.L2.AND.NTYP.EQ.3) THEN
C
C     THREE POINT QUADRATIC EXTRAPOLATION OR INTERPOLATION
C
      AL=(FI-FL)/(RI-RL)
      AU=(FU-FI)/(RU-RI)
      DEN=AU-AL
      IF (DEN.LE.ZERO) GO TO 1
      R=RI-HALF*(AU*(RI-RL)+AL*(RU-RI))/DEN
      ELSE IF (L1.AND..NOT.L2.AND.NTYP.EQ.4) THEN
C
C     THREE POINT CUBIC EXTRAPOLATION OR INTERPOLATION
C
      DIS=(AI-ONE)*(RU/RI)
      DEN=(AU-ONE)*(RI/RU)-DIS
      DIS=AU+AI-DEN-TWO*(ONE+DIS)
      DIS=DEN*DEN-THREE*DIS
      IF (DIS.LT.ZERO) GO TO 1
      DEN=DEN+SQRT(DIS)
      IF (DEN.EQ.ZERO) GO TO 1
      R=(RU-RI)/DEN
      ELSE
      GO TO 1
      END IF
      IF (MODE .EQ. 1  .AND.  R .GT. RU)  THEN
C
C     EXTRAPOLATION ACCEPTED
C
      R = MAX( R, C1L*RU)
      R = MIN( R, C1U*RU)
      RETURN
      ELSE IF (MODE .EQ. 2 .AND. R .GT. RL .AND. R .LT. RU) THEN
C
C     INTERPOLATION ACCEPTED
C
      IF (RI.EQ.ZERO.AND.NTYP.NE.4) THEN
      R = MAX( R, RL + C2L*(RU-RL))
      ELSE
      R = MAX( R, RL + C3L*(RU-RL))
      END IF
      R = MIN( R, RL + C2U*(RU-RL))
      IF (R.EQ.RI) GO TO 1
      RETURN
      END IF
    1 CONTINUE
      END
* SUBROUTINE PP0AF1             ALL SYSTEMS                 97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* COMPUTATION OF VALUE OF THE AUGMENTED LAGRANGIAN FUNCTION.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  II  NC  NUMBER OF CONSTRAINTS.
*  RI  CF(NC+1)  VECTOR CONTAINING VALUES OF THE CONSTRAINTS.
*  RI  CZ(NC)  VECTOR OF LAGRANGE MULTIPLIERS.
*  RI  S(NF+NC)  DIRECTION VECTOR.
*  RI  RPF2  PENALTY COEFFICIENT.
*  RO  FC  VALUE OF THE PENALTY TERM.
*  RO  F  VALUE OF THE PENALTY FUNCTION.
*
      SUBROUTINE PP0AF1(NF,NC,CF,CZ,S,RPF2,FC,F)
      INTEGER NF,NC
      DOUBLE PRECISION CF(NC+1),CZ(NC),S(NF+NC),RPF2,FC,F
      INTEGER KC
      FC=0.0D 0
      DO 1 KC=1,NC
      FC=FC+(CZ(KC)+S(NF+KC)+0.5D 0*RPF2*CF(KC))*CF(KC)
    1 CONTINUE
      F=CF(NC+1)+FC
      RETURN
      END
* SUBROUTINE PP1LF3             ALL SYSTEMS                 97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* COMPUTATION OF THE VALUE AND THE GRADIENT OF THE OBJECTIVE FUNCTION.
* COMPUTATION OF THE VALUES AND THE GRADIENTS OF THE CONSTRAINT
* FUNCTIONS.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  II  NC  NUMBER OF CONSTRAINTS.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  RI  GF(NF)  GRADIENT OF THE MODEL FUNCTION.
*  RO  CF(NC)  VECTOR OF THE CONSTRAINT FUNCTION VALUES.
*  RI  GC(NF)  GRADIENT OF THE SELECTED FUNCTION.
*  II  ICG(NC+1)  POSITIONS OF THE FIRST ROWS ELEMENTS IN THE SPARSE
*         JACOBIAN STRUCTURE.
*  II  JCG(MC)  COLUMN INDICES OF ELEMENTS IN THE SPARSE JACOBIAN
*         STRUCTURE.
*  RA  GC(NF) GRADIENT OF THE CONSTRAINT FUNCTION.
*  RO  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RI  CZ(NC)  VECTOR OF LAGRANGE MULTIPLIERS.
*  RI  FF  VALUE OF THE MODEL FUNCTION.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RO  CMAX  MAXIMUM CONSTRAINT VIOLATION.
*  IU  KD  DEGREE OF REQUIRED DERVATIVES.
*  IU  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES.
*  IU  NFV  NUMBER OF OBJECTIVE FUNCTION VALUES COMPUTED.
*  IU  NFG  NUMBER OF OBJECTIVE FUNCTION GRADIENTS COMPUTED.
*  II  ISNC  INDICATOR FOR STORING CONSTRAINT FUNCTION VALUES AND
*         GRADIENTS. ISNC=0-STORING SUPPRESSED. ISNC=1-STORING
*         CONSTRAINT FUNCTION VALUES. ISNC=2-STORING CONSTRAINT
*         FUNCTION VALUES AND GRADIENTS.
*
* SUBPROGRAMS USED :
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVNEG  COPYING OF A VECTOR WITH CHANGE OF THE SIGN.
*
      SUBROUTINE PP1LF3(NF,NC,X,GF,CF,CG,ICG,JCG,GC,G,CZ,FF,F,CMAX,
     & KD,LD,NFV,NFG,ISNC)
      INTEGER NF,NC,ICG(*),JCG(*),KD,LD,NFV,NFG,ISNC
      DOUBLE PRECISION X(*),GF(*),CF(*),CG(*),GC(*),G(*),CZ(*),FF,F,CMAX
      INTEGER J,K,L,JP,KC
      DOUBLE PRECISION FC
      IF (KD.LE.LD) RETURN
      IF (LD.GE.0) GO TO 1
      CALL OBJ(NF,X,FF)
      NFV=NFV+1
      F=FF
      CMAX=0.0D 0
    1 IF (KD.LT.1) GO TO 2
      IF (LD.GE.1) GO TO 2
      CALL DOBJ(NF,X,GF)
      NFG=NFG+1
      CALL MXVCOP(NF,GF,G)
    2 DO 5 KC=1,NC
      IF (LD.GE.0) GO TO 3
      CALL CON(NF,KC,X,FC)
      IF (ISNC.GE.0) CF(KC)=FC
      F=F+CZ(KC)*FC
      CMAX=MAX(CMAX,ABS(FC))
    3 IF (KD.LT.1) GO TO 5
      IF (LD.GE.1) GO TO 5
      CALL DCON(NF,KC,X,GC)
      K=ICG(KC)
      L=ICG(KC+1)-K
      DO 4 J=1,L
      JP=JCG(K)
      G(JP)=G(JP)+CZ(KC)*GC(JP)
      IF (ISNC.GE.1) CG(K)=GC(JP)
      K=K+1
    4 CONTINUE
    5 CONTINUE
      LD=KD
      RETURN
      END
* SUBROUTINE PS0L02                ALL SYSTEMS                97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
*  EXTENDED LINE SEARCH WITHOUT DIRECTIONAL DERIVATIVES.
*
* PARAMETERS :
*  RO  R  VALUE OF THE STEPSIZE PARAMETER.
*  RO  RO  INITIAL VALUE OF THE STEPSIZE PARAMETER.
*  RO  RP  PREVIOUS VALUE OF THE STEPSIZE PARAMETER.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  INITIAL VALUE OF THE OBJECTIVE FUNCTION.
*  RO  FP  PREVIOUS VALUE OF THE OBJECTIVE FUNCTION.
*  RI  PO  INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RO  PP  PREVIOUS VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RI  FMIN  LOWER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FMAX  UPPER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION.
*  RI  RMIN  MINIMUM VALUE OF THE STEPSIZE PARAMETER
*  RI  RMAX  MAXIMUM VALUE OF THE STEPSIZE PARAMETER
*  RI  GNORM  NORM OF THE GRADIENT VECTOR.
*  RI  SNORM  NORM OF THE DIRECTION VECTOR.
*  RI  TOLS  TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE
*         CHANGE OF THE FUNCTION VALUE).
*  II  NIT  ACTUAL NUMBER OF ITERATIONS.
*  II  KIT  NUMBER OF THE ITERATION AFTER LAST RESTART.
*  II  IEST  LOWER BOUND SPECIFICATION. IEST=0 OR IEST=1 IF LOWER BOUND
*         IS NOT OR IS GIVEN.
*  II  INITS  CHOICE OF THE INITIAL STEPSIZE. INITS=0-INITIAL STEPSIZE
*         IS SPECIFIED IN THE CALLING PROGRAM. INITS=1-UNIT INITIAL
*         STEPSIZE. INITS=2-COMBINED UNIT AND QUADRATICALLY ESTIMATED
*         INITIAL STEPSIZE. INITS=3-QUADRATICALLY ESTIMATED INITIAL
*         STEPSIZE.
*  IO  ITERS  TERMINATION INDICATOR. ITERS=0-ZERO STEP. ITERS=1-PERFECT
*         LINE SEARCH. ITERS=2 GOLDSTEIN STEPSIZE. ITERS=3-CURRY
*         STEPSIZE. ITERS=4-EXTENDED CURRY STEPSIZE.
*         ITERS=5-ARMIJO STEPSIZE. ITERS=6-FIRST STEPSIZE.
*         ITERS=7-MAXIMUM STEPSIZE. ITERS=8-UNBOUNDED FUNCTION.
*         ITERS=-1-MRED REACHED. ITERS=-2-POSITIVE DIRECTIONAL
*         DERIVATIVE. ITERS=-3-ERROR IN INTERPOLATION.
*  IO  MAXST  MAXIMUM STEPSIZE INDICATOR. MAXST=0 OR MAXST=1 IF MAXIMUM
*         STEPSIZE WAS NOT OR WAS REACHED.
*  IO  NRED  ACTUAL NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS.
*
*  II  MRED  MAXIMUM NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS.
*  II  KTERS  TERMINATION SELECTION. KTERS=1-PERFECT LINE SEARCH.
*         KTERS=2-GOLDSTEIN STEPSIZE. KTERS=3-CURRY STEPSIZE.
*         KTERS=4-EXTENDED CURRY STEPSIZE. KTERS=5-ARMIJO STEPSIZE.
*         KTERS=6-FIRST STEPSIZE.
*  II  MES  METHOD SELECTION. MES=1-BISECTION. MES=2-QUADRATIC
*         INTERPOLATION (WITH ONE DIRECTIONAL DERIVATIVE).
*         MES=3-QUADRATIC INTERPOLATION (WITH TWO DIRECTIONAL
*         DERIVATIVES). MES=4-CUBIC INTERPOLATION. MES=5-CONIC
*         INTERPOLATION.
*  II  MES1  SWITCH FOR EXTRAPOLATION. MES1=1-CONSTANT INCREASING OF
*         THE INTERVAL. MES1=2-EXTRAPOLATION SPECIFIED BY THE PARAMETER
*         MES. MES1=3 SUPPRESSED EXTRAPOLATION.
*  II  MES2  SWITCH FOR TERMINATION. MES2=1-NORMAL TERMINATION.
*         MES2=2-TERMINATION AFTER AT LEAST TWO STEPS (ASYMPTOTICALLY
*         PERFECT LINE SEARCH).
*  II  MES3  SAFEGUARD AGAINST ROUNDING ERRORS. MES3=0-SAFEGUARD
*         SUPPRESSED. MES3=1-FIRST LEVEL OF SAFEGUARD. MES3=2-SECOND
*         LEVEL OF SAFEGUARD.
*  IU  ISYS  CONTROL PARAMETER.
*
* SUBPROGRAM USED :
*  S   PNINT3  EXTRAPOLATION OR INTERPOLATION WITHOUT DIRECTIONAL
*         DERIVATIVES.
*
* METHOD :
* SAFEGUARDED EXTRAPOLATION AND INTERPOLATION WITH EXTENDED TERMINATION
* CRITERIA.
*
      SUBROUTINE PS0L02(R,RO,RP,F,FO,FP,PO,PP,FMIN,FMAX,RMIN,RMAX,
     & TOLS,KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,INITS,ITERS,KTERS,
     & MES,ISYS)
      INTEGER KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,INITS,ITERS,KTERS,
     & MES,ISYS
      DOUBLE PRECISION R,RO,RP,F,FO,FP,PO,PP,FMIN,FMAX,RMIN,RMAX,TOLS
      DOUBLE PRECISION RL,FL,RU,FU,RI,FI,RTEMP,TOL
      INTEGER MTYP,MERR,MODE,INIT1,MES1,MES2
      LOGICAL L1,L2,L3,L4,L6,L7
      PARAMETER(TOL=1.0D-4)
      SAVE MTYP,MODE,MES1,MES2
      SAVE RL,FL,RU,FU,RI,FI
      IF (ISYS.EQ.1) GO TO 3
C      GO TO (1,3) ISYS+1
    1 CONTINUE
      MES1=2
      MES2=2
      ITERS=0
      IF (PO.GE.0.0D 0) THEN
      R=0.0D 0
      ITERS=-2
      GO TO 4
      END IF
      IF(RMAX.LE.0.0D 0) THEN
      ITERS= 0
      GO TO 4
      END IF
C
C     INITIAL STEPSIZE SELECTION
C
      IF (INITS.GT.0) THEN
      RTEMP=FMIN-F
      ELSE IF (IEST.EQ.0) THEN
      RTEMP=F-FP
      ELSE
      RTEMP=MAX(F-FP,1.0D 1*(FMIN-F))
      END IF
      INIT1=ABS(INITS)
      RP=0.0D 0
      FP=FO
      PP=PO
      IF (INIT1.EQ.0) THEN
      ELSE IF (INIT1.EQ.1.OR.INITS.GE.1.AND.IEST.EQ.0) THEN
      R=1.0D 0
      ELSE IF (INIT1.EQ.2) THEN
      R=MIN(1.0D 0,4.0D 0*RTEMP/PO)
      ELSE IF (INIT1.EQ.3) THEN
      R=MIN(1.0D 0, 2.0D 0*RTEMP/PO)
      ELSE IF (INIT1.EQ.4) THEN
      R=2.0D 0*RTEMP/PO
      END IF
      RTEMP=R
      R=MAX(R,RMIN)
      R=MIN(R,RMAX)
      MODE=0
      RL=0.0D 0
      FL=FO
      RU=0.0D 0
      FU=FO
      RI=0.0D 0
      FI=FO
C
C     NEW STEPSIZE SELECTION (EXTRAPOLATION OR INTERPOLATION)
C
    2 CALL PNINT3(RO,RL,RU,RI,FO,FL,FU,FI,PO,R,MODE,MTYP,MERR)
      IF (MERR.GT.0) THEN
      ITERS=-MERR
      GO TO 4
      ELSE IF (MODE.EQ.1) THEN
      NRED=NRED-1
      R=MIN(R,RMAX)
      ELSE IF (MODE.EQ.2) THEN
      NRED=NRED+1
      END IF
C
C     COMPUTATION OF THE NEW FUNCTION VALUE
C
      KD= 0
      LD=-1
      ISYS=1
      RETURN
    3 CONTINUE
      IF (ITERS.NE.0) GO TO 4
      IF (F.LE.FMIN) THEN
      ITERS=7
      GO TO 4
      ELSE
      L1=R.LE.RMIN.AND.NIT.NE.KIT
      L2=R.GE.RMAX
      L3=F-FO.LE.TOLS*R*PO.OR.F-FMIN.LE.(FO-FMIN)/1.0D 1
      L4=F-FO.GE.(1.0D 0-TOLS)*R*PO.OR.MES2.EQ.2.AND.MODE.EQ.2
      L6=RU-RL.LE.TOL*RU.AND.MODE.EQ.2
      L7=MES2.LE.2.OR.MODE.NE.0
      MAXST=0
      IF (L2) MAXST=1
      END IF
C
C     TEST ON TERMINATION
C
      IF (L1.AND..NOT.L3) THEN
      ITERS=0
      GO TO 4
      ELSE IF (L2.AND..NOT.F.GE.FU) THEN
      ITERS=7
      GO TO 4
      ELSE IF (L6) THEN
      ITERS=1
      GO TO 4
      ELSE IF (L3.AND.L7.AND.KTERS.EQ.5) THEN
      ITERS=5
      GO TO 4
      ELSE IF (L3.AND.L4.AND.L7.AND.(KTERS.EQ.2.OR.KTERS.EQ.3.OR.
     *         KTERS.EQ.4)) THEN
      ITERS=2
      GO TO 4
      ELSE IF (KTERS.LT.0.OR.KTERS.EQ.6.AND.L7) THEN
      ITERS=6
      GO TO 4
      ELSE IF(ABS(NRED).GE.MRED) THEN
      ITERS=-1
      GO TO 4
      ELSE
      RP=R
      FP=F
      MODE=MAX(MODE,1)
      MTYP=ABS(MES)
      IF(F.GE.FMAX) MTYP=1
      END IF
      IF (MODE.EQ.1) THEN
C
C     INTERVAL CHANGE AFTER EXTRAPOLATION
C
      RL=RI
      FL=FI
      RI=RU
      FI=FU
      RU=R
      FU=F
      IF (F.GE.FI) THEN
      NRED=0
      MODE=2
      ELSE IF ( MES1 .EQ. 1) THEN
      MTYP=1
      END IF
C
C     INTERVAL CHANGE AFTER INTERPOLATION
C
      ELSE IF (R.LE.RI) THEN
      IF (F.LE.FI) THEN
      RU=RI
      FU=FI
      RI=R
      FI=F
      ELSE
      RL=R
      FL=F
      END IF
      ELSE
      IF (F.LE.FI) THEN
      RL=RI
      FL=FI
      RI=R
      FI=F
      ELSE
      RU=R
      FU=F
      END IF
      END IF
      GO TO 2
    4 ISYS=0
      RETURN
      END
* SUBROUTINE PYCSER                ALL SYSTEMS                98/12/01
C PORTABILITY : ALL SYSTEMS
C 98/12/01 TU : ORIGINAL VERSION
*
* PURPOSE :
* GROUP OF THE SAME COLOUR FOR THE POWELL-TOINT ALGORITHM FOR SPARSE
* HESSIANS APPROXIMATIONS IS CREATED.
*
* PARAMETERS :
*  IU  IH(MCOLS+1) POINTER VECTOR OF SPARSE HESSIAN MATRIX.
*  IU  JH(M) INDEX VECTOR OF THE HESSIAN MATRIX.
*  IA  WN02(MCOLS) AUXILIARY VECTOR.
*  RA  WN03(MCOLS) AUXILIARY VECTOR.
*  RI  DEG(MCOLS) DEGREES OF THE ADJACENCY GRAPH.
*  IA  WN01(NF) AUXILIARY VECTOR USED FOR INDICES OF THE COLUMNS
*        THAT HAVE NOT BEEN COLOURED YET.
*  II  COL(NF) VECTOR DISCERNING GROUPS OF THE HESSIAN COLUMN OF THE
*              SAME COLOUR.
*  IU  NCOL  NUMBER OF COLOURS USED SO FAR.
*  IU  CNM  NUMBER OF COLUMNS THAT HAVE NOT BEEN COLOURED SO FAR.
*
      SUBROUTINE PYCSER(JH,IH,WN02,WN03,DEG,WN01,COL,NCOL,CNM)
      INTEGER JH(*),IH(*),COL(*)
      INTEGER WN01(*),WN02(*)
      DOUBLE PRECISION WN03(*),DEG(*)
      INTEGER NCOL,CNM,I,J,K,L,IP
C
C     DEFINITION OF THE INCIDENCE ARRAY A
C
      L=WN01(1)
C
C     ELEMENT WAS MARKED THAT IT IS INSERTED
C
      DO 100 I=IH(L),IH(L+1)-1
      K=JH(I)
C
C     COLUMN OF THIS NUMBER HAS APPEARED IN ONE OF THE PREVIOUS GROUPS
C
      IF(COL(K).LT.NCOL) GO TO 100
      DEG(K)=DEG(K)-1
      WN02(K)=NCOL
100   CONTINUE
C
C     COLUMN IS INSERTED
C
      COL(L)=NCOL
C
C     THE CYCLE OF COMPARING COLUMN WITH THE ARRAY A
C     A2 IS AN HELP ARRAY CONTAINING COLUMNS THAT ARE
C     BEEING EXAMINED BUT THAT WERE NOT YET ACCEPTED
C     P IS ITS POINTER
C
      IF(CNM.EQ.1) GO TO 250
      DO 200 I=2,CNM
C
C     TRANSFORMATION OF THE EXAMINED COLUMN I IS
C
      IP=1
      L=WN01(I)
      DO 300 J=IH(L),IH(L+1)-1
        K=JH(J)
        IF(COL(K).LT.NCOL) GO TO 300
        IF(WN02(K).GE.NCOL) GO TO 200
C
        WN03(IP)=K
        IP=IP+1
300   CONTINUE
      IF(IP.NE.1) THEN
C
C COPY OF THE WN03 ARRAY INTO WN02 FOR THE COLUMN WAS ACCEPTED
C
      DO 400 K=1,IP-1
        WN02(INT(WN03(K)))=NCOL
        DEG(INT(WN03(K)))=DEG(INT(WN03(K)))-1
400   CONTINUE
      END IF
C
C     INSERT THE COLUMN INTO THE PROCESSED GROUP
C
      COL(L)=NCOL
C
C     END OF THE MAIN CYCLE
C
200   CONTINUE
C
C     JUMP LABEL
C
250   CONTINUE
C
C     INVP SHIFT
C
      K=1
      DO 500 I=1,CNM
        L=WN01(I)
        IF(COL(L).EQ.NCOL) THEN
      ELSE
        WN01(K)=L
        K=K+1
      END IF
500   CONTINUE
C
C     CNM UPDATE
C
      CNM=K-1
C
      RETURN
      END
* SUBROUTINE PYPTSH                ALL SYSTEMS                98/12/01
C PORTABILITY : ALL SYSTEMS
C 98/12/01 TU : ORIGINAL VERSION
*
* PURPOSE :
* POWELL-TOINT GRAPH COLORING ALGORITHM FOR GROUPING COLUMNS OF THE
* HESSIAN MATRIX BEFORE NUMERICAL DIFFERENTIATION.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  IH(NF+1) POINTER VECTOR OF SPARSE HESSIAN MATRIX.
*  II  JH(MMAX) INDEX VECTOR OF THE HESSIAN MATRIX.
*  IO  COL(NF) VECTOR DISCERNING GROUPS OF THE HESSIAN COLUMN OF THE
*              SAME COLOUR.
*  RA  DEG(NF) DEGREES OF THE ADJACENCY GRAPH.
*  RA  ORD(NF) AUXILIARY ARRAY.
*  RA  RADIX(NF+1) AUXILIARY ARRAY.
*  IA  WN11(NF) AUXILIARY VECTOR USED FOR INDICES OF THE COLUMNS
*        THAT HAVE NOT BEEN COLOURED YET.
*  IA  WN12(NF) AUXILIARY VECTOR.
*  RA  XS(NF) AUXILIARY VECTOR.
*
* SUBPROGRAMS USED :
*  S   PYCSER  GROUPING COLUMNS OF THE SPARSE SYMMETRIC MATRIX.
*  S   MXSTG1  WIDTHEN THE STRUCTURE.
*  S   MXSTL1  SHRINK THE STRUCTURE.
*  S   MXVSR2  SORT.
*
      SUBROUTINE PYPTSH(NF,MMAX,IH,JH,COL,DEG,ORD,RADIX,WN11,WN12,XS,
     & ITERM)
      INTEGER NF,MMAX,IH(*),JH(*),COL(*)
      INTEGER WN11(*),WN12(*),ITERM
      DOUBLE PRECISION RADIX(*),ORD(*)
      DOUBLE PRECISION XS(*),DEG(*)
      INTEGER NCOL,CNM,I,ML,MM,J,K1,L
C
C SAVE SYMBOLIC STRUCTURE OF FACTOR
C
      MM=IH(NF+1)-1
      IF(2*MM-NF+2.GE.MMAX) THEN
        ITERM=-45
        RETURN
      END IF
C
C WIDTHEN THE STRUCTURE
C
      CALL MXSTG1(NF,ML,IH,JH,WN12,WN11)
C
      DO 100 I=1,NF
      COL(I)=NF
      WN12(I)=0
      WN11(I)=I
100   CONTINUE
C
C     NUMBER OF THE FREE COLUMNS
C
      CNM=NF
C
C     NUMBER OF USED COLOURS
C
      NCOL=1
C
C     DEGREE RECOUNT
C
      K1=1
      DO 110 I=1,NF
      L=IH(I+1)
      DEG(I)=L-K1
      K1=L
110   CONTINUE
C
C     COLUMN RESORT
C
200   CALL MXVSR2(NF,DEG,ORD,RADIX,WN11,CNM)
C
C     ORD REWRITE INTO THE ARRAY INVP
C
      DO 250 I=1,CNM
        WN11(I)=ORD(I)
250   CONTINUE
C
C     COLUMNS OF THE NEW COLOUR NCOL
C
      CALL PYCSER(JH,IH,WN12,XS,DEG,WN11,COL,NCOL,CNM)
C
C     STOP TEST
C
      IF(CNM.GE.1) THEN
        NCOL=NCOL+1
        GO TO 200
      END IF
C
C SHRINK THE STRUCTURE
C
      CALL MXSTL1(NF,ML,IH,JH,WN12)
C
C INTO COL GIVE INDICES OF THE INDIVIDUAL GROUPS ONE AFTER ANOTHER,
C END OF THE GROUP IS MARKED BY THE NEGATIVE INDEX VALUE.
C
C
C READ COL
C
      DO 300 I=1,NF
        WN11(I)=0
 300  CONTINUE
      DO 400 I=1,NF
        J=COL(I)
        WN11(J)=WN11(J)+1
 400  CONTINUE
      WN12(1)=1
      L=1
      DO 500 I=2,NF
        L=L+WN11(I-1)
        WN12(I)=L
        IF(WN11(I).EQ.0) GO TO 550
 500  CONTINUE
 550  CONTINUE
C
C CHANGE COL
C
      DO 600 I=1,NF
        J=COL(I)
        WN11(I)=J
 600  CONTINUE
      DO 700 I=1,NF
        J=WN11(I)
        COL(WN12(J))=I
        WN12(J)=WN12(J)+1
 700  CONTINUE
      DO 800 I=1,NCOL
        L=WN12(I)-1
        IF(L.GT.NF) GO TO 900
        COL(L)=-COL(L)
 800  CONTINUE
 900  CONTINUE
      RETURN
      END