************************************************************************
* SUBROUTINE PNEDU              ALL SYSTEMS                   97/01/22
* PURPOSE :
* EASY TO USE SUBROUTINE FOR UNCONSTRAINED MINIMIZATION OF
* FUNCTIONS WITH LARGE-SCALE SPARSE HESSIAN MATRICES
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  IU  MH  NUMBER OF NONZERO ELEMENTS IN THE UPPER PART OF THE HESSIAN
*         MATRIX.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  II  IH(NF+1)  POINTERS OF DIAGONAL ELEMENTS OF THE MATRIX H.
*  IU  JH(MH)  INDICES OF NONZERO ELEMENTS OF THE MATRIX H
*         TOGETHER WITH AN ADDITIONAL SPACE USED FOR THE NUMERICAL
*         DIFFERENTIATION.
*  II  IPAR(7)  INTEGER PAREMETERS:
*      IPAR(1)  MAXIMUM NUMBER OF ITERATIONS.
*      IPAR(2)  MAXIMUM NUMBER OF FUNCTION EVALUATIONS.
*      IPAR(3)  MAXIMUM NUMBER OF GRADIENT EVALUATIONS.
*      IPAR(4)  ESTIMATION INDICATOR. IPAR(4)=0-MINIMUM IS NOT
*         ESTIMATED. IPAR(4)=1-MINIMUM IS ESTIMATED BY THE VALUE
*         RPAR(6).
*      IPAR(5)  METHOD FOR COMPUTING THE TRUST REGION STEP.
*         IPAR(5)=1-THE DOG-LEG METHOD. MOS=2-THE MORE-SORENSEN METHOD.
*      IPAR(6)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PNED.
*      IPAR(7)  NUMBER DEFINING THE SPACE FOR FILL-IN (THE SIZE OF THIS
*         SPACE IS IFIL TIMES THE STANDARD SIZE). THE DEFAULT VALUE IS
*         IFIL=1. THE DEFAULT VALUE HAS TO BE INCREASED IF ITERM IS
*         LESS OR EQUAL TO -40.
*  RI  RPAR(9)  REAL PARAMETERS:
*      RPAR(1)  MAXIMUM STEPSIZE.
*      RPAR(2)  TOLERANCE FOR THE CHANGE OF VARIABLES.
*      RPAR(3)  TOLERANCE FOR THE CHANGE OF FUNCTION VALUES.
*      RPAR(4)  TOLERANCE FOR THE FUNCTION FALUE.
*      RPAR(5)  TOLERANCE FOR THE GRADIENT NORM.
*      RPAR(6)  ESTIMATION OF THE MINIMUM FUNCTION VALUE.
*      RPAR(7)  INITIAL TRUST-REGION RADIUS.
*      RPAR(8)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PNED.
*      RPAR(9)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PNED.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RO  GMAX  MAXIMUM PARTIAL DERIVATIVE.
*  II  ISPAS  INPUT SPARSE STRUCTURE. ISPAS=1-STANDARD COORDINATE
*         FORM. ISPAS=2-SPARSE STRUCTURE COMPRESSED BY ROWS.
*  II  IPRNT  PRINT SPECIFICATION. IPRNT=0-NO PRINT.
*         ABS(IPRNT)=1-PRINT OF FINAL RESULTS.
*         ABS(IPRNT)=2-PRINT OF FINAL RESULTS AND ITERATIONS.
*         IPRNT>0-BASIC FINAL RESULTS. IPRNT<0-EXTENDED FINAL
*         RESULTS.
*  IO  ITERM  VARIABLE THAT INDICATES THE CAUSE OF TERMINATION.
*         ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN
*                   MTESX (USUALLY TWO) SUBSEQUENT ITERATIONS.
*         ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN
*                   MTESF (USUALLY TWO) SUBSEQUENT ITERATIONS.
*         ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB.
*         ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG.
*         ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED,
*                   BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE.
*         ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV.
*         ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED.
*         VALUES ITERM<=-40 DETECT A LACK OF SPACE. IN THIS CASE,
*         PARAMETER IPAR(7) HAS TO BE INCREASED.
*
* VARIABLES IN COMMON /STAT/ (STATISTICS) :
*  IO  NRES  NUMBER OF RESTARTS.
*  IO  NDEC  NUMBER OF MATRIX DECOMPOSITIONS.
*  IO  NIN  NUMBER OF INNER ITERATIONS.
*  IO  NIT  NUMBER OF ITERATIONS.
*  IO  NFV  NUMBER OF FUNCTION EVALUATIONS.
*  IO  NFG  NUMBER OF GRADIENT EVALUATIONS.
*  IO  NFH  NUMBER OF HESSIAN EVALUATIONS.
*
* SUBPROGRAMS USED :
*  S   PNED  DISCRETE NEWTON METHOD WITH DIRECT SOLUTION OF THE
*         TRUST-REGION SUBPROBLEM.
*  S   PFSED3  COMPRESSED SPARSE STRUCTURE OF THE HESSIAN MATRIX IS
*         COMPUTED FROM THE COORDINATE FORM.
*  S   MXVICP COPYING OF AN INTEGER VECTOR.
*
* EXTERNAL SUBROUTINES :
*  SE  OBJ  COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION.
*         CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS THE NUMBER
*         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND FF IS
*         THE VALUE OF THE OBJECTIVE FUNCTION.
*  SE  DOBJ  COMPUTATION OF THE GRADIENT OF THE OBJECTIVE FUNCTION.
*         CALLING SEQUENCE: CALL DOBJ(NF,X,GF) WHERE NF IS THE NUMBER
*         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND GF(NF)
*         IS THE GRADIENT OF THE OBJECTIVE FUNCTION.
*
      SUBROUTINE PNEDU(NF,MH,X,IH,JH,IPAR,RPAR,F,GMAX,ISPAS,IPRNT,
     & ITERM)
      INTEGER NF,MH,IH(*),JH(*),IPAR(7),ISPAS,IPRNT,ITERM
      DOUBLE PRECISION X(*),RPAR(9),F,GMAX
      INTEGER NB,LGF,LHF,LS,LXO,LGO,LXS,LCOL,LPSL,LPERM,LINVP,LWN11,
     & LWN12,LWN13,LWN14,LJH,IFIL,IER
      INTEGER NRES,NDEC,NIN,NIT,NFV,NFG,NFH
      COMMON /STAT/ NRES,NDEC,NIN,NIT,NFV,NFG,NFH
      INTEGER IA(:)
      DOUBLE PRECISION RA(:)
      ALLOCATABLE IA,RA
      IFIL=IPAR(7)
      IF (IFIL.LE.0) IFIL=1
      IF (ISPAS.LE.1) THEN
      CALL PFSED3(NF,MH,IH,JH,IER)
      IF (IER.NE.0) THEN
      WRITE (6,'(''INPUT ERROR : IER = '',I3)') IER
      STOP
      END IF
      ELSE
      MH=IH(NF+1)-1
      END IF
      ALLOCATE (IA(8*NF+5+(2*IFIL+3)*MH),RA(5*NF+1+(2*IFIL+3)*MH))
      NB=0
*
*     POINTERS FOR AUXILIARY ARRAYS
*
      LGF=1
      LS=LGF+NF
      LXO=LS+NF
      LGO=LXO+NF
      LXS=LGO+NF+1
      LHF=LXS+NF
      LCOL=1
      LPSL=LCOL+NF
      LPERM=LPSL+NF+1
      LINVP=LPERM+NF
      LWN11=LINVP+NF
      LWN12=LWN11+NF+1
      LWN13=LWN12+NF+1
      LWN14=LWN13+NF+1
      LJH=LWN14+NF+1
      CALL MXVICP(IH(NF+1)-1,JH,IA(LJH))
      CALL PNED(NF,NB,(2*IFIL+3)*MH,X,IA,RA,RA,RA(LGF),RA(LHF),IH,
     & IA(LJH),RA(LS),RA(LXO),RA(LGO),RA(LXS),IA(LCOL),IA(LPSL),
     & IA(LPERM),IA(LINVP),IA(LWN11),IA(LWN12),IA(LWN13),IA(LWN14),
     & RPAR(1),RPAR(2),RPAR(3),RPAR(4),RPAR(5),RPAR(6),RPAR(7),GMAX,
     & F,IPAR(1),IPAR(2),IPAR(3),IPAR(4),IPAR(5),IPRNT,ITERM)
      DEALLOCATE (IA,RA)
      RETURN
      END
************************************************************************
* SUBROUTINE PNEDS              ALL SYSTEMS                   97/01/22
* PURPOSE :
* EASY TO USE SUBROUTINE FOR BOX CONSTRAINED MINIMIZATION OF
* LARGE-SCALE FUNCTIONS WITH SPARSE HESSIAN MATRICES
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  IU  MH  NUMBER OF NONZERO ELEMENTS IN THE UPPER PART OF THE HESSIAN
*         MATRIX.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS. IX(I)=0-VARIABLE
*         X(I) IS UNBOUNDED. IX(I)=1-LOWER BOUND XL(I).LE.X(I).
*         IX(I)=2-UPPER BOUND X(I).LE.XU(I). IX(I)=3-TWO SIDE BOUND
*         XL(I).LE.X(I).LE.XU(I). IX(I)=5-VARIABLE X(I) IS FIXED.
*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
*  II  IH(NF+1)  POINTERS OF DIAGONAL ELEMENTS OF THE MATRIX H.
*  IU  JH(MH)  INDICES OF NONZERO ELEMENTS OF THE MATRIX H
*         TOGETHER WITH AN ADDITIONAL SPACE USED FOR THE NUMERICAL
*         DIFFERENTIATION.
*  II  IPAR(7)  INTEGER PAREMETERS:
*      IPAR(1)  MAXIMUM NUMBER OF ITERATIONS.
*      IPAR(2)  MAXIMUM NUMBER OF FUNCTION EVALUATIONS.
*      IPAR(3)  MAXIMUM NUMBER OF GRADIENT EVALUATIONS.
*      IPAR(4)  ESTIMATION INDICATOR. IPAR(4)=0-MINIMUM IS NOT
*         ESTIMATED. IPAR(4)=1-MINIMUM IS ESTIMATED BY THE VALUE
*         RPAR(6).
*      IPAR(5)  METHOD FOR COMPUTING THE TRUST REGION STEP. MOS=1-THE
*         DOG-LEG METHOD. MOS=2-THE MORE-SORENSEN METHOD.
*      IPAR(6)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PNED.
*      IPAR(7)  NUMBER DEFINING THE SPACE FOR FILL-IN (THE SIZE OF THIS
*         SPACE IS IFIL TIMES THE STANDARD SIZE). THE DEFAULT VALUE IS
*         IFIL=1. THE DEFAULT VALUE HAS TO BE INCREASED IF ITERM IS
*         LESS OR EQUAL TO -40.
*  RI  RPAR(9)  REAL PARAMETERS:
*      RPAR(1)  MAXIMUM STEPSIZE.
*      RPAR(2)  TOLERANCE FOR THE CHANGE OF VARIABLES.
*      RPAR(3)  TOLERANCE FOR THE CHANGE OF FUNCTION VALUES.
*      RPAR(4)  TOLERANCE FOR THE FUNCTION FALUE.
*      RPAR(5)  TOLERANCE FOR THE GRADIENT NORM.
*      RPAR(6)  ESTIMATION OF THE MINIMUM FUNCTION VALUE.
*      RPAR(7)  INITIAL TRUST-REGION RADIUS.
*      RPAR(8)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PNED.
*      RPAR(9)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PNED.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RO  GMAX  MAXIMUM PARTIAL DERIVATIVE.
*  II  ISPAS  INPUT SPARSE STRUCTURE. ISPAS=1-STANDARD COORDINATE
*         FORM. ISPAS=2-SPARSE STRUCTURE COMPRESSED BY ROWS.
*  II  IPRNT  PRINT SPECIFICATION. IPRNT=0-NO PRINT.
*         ABS(IPRNT)=1-PRINT OF FINAL RESULTS.
*         ABS(IPRNT)=2-PRINT OF FINAL RESULTS AND ITERATIONS.
*         IPRNT>0-BASIC FINAL RESULTS. IPRNT<0-EXTENDED FINAL
*         RESULTS.
*  IO  ITERM  VARIABLE THAT INDICATES THE CAUSE OF TERMINATION.
*         ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN
*                   MTESX (USUALLY TWO) SUBSEQUENT ITERATIONS.
*         ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN
*                   MTESF (USUALLY TWO) SUBSEQUENT ITERATIONS.
*         ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB.
*         ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG.
*         ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED,
*                   BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE.
*         ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV.
*         ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED.
*         VALUES ITERM<=-40 DETECT A LACK OF SPACE. IN THIS CASE,
*         PARAMETER IPAR(7) HAS TO BE INCREASED.
*
* VARIABLES IN COMMON /STAT/ (STATISTICS) :
*  IO  NRES  NUMBER OF RESTARTS.
*  IO  NDEC  NUMBER OF MATRIX DECOMPOSITIONS.
*  IO  NIN  NUMBER OF INNER ITERATIONS.
*  IO  NIT  NUMBER OF ITERATIONS.
*  IO  NFV  NUMBER OF FUNCTION EVALUATIONS.
*  IO  NFG  NUMBER OF GRADIENT EVALUATIONS.
*  IO  NFH  NUMBER OF HESSIAN EVALUATIONS.
*
* SUBPROGRAMS USED :
*  S   PNED  DISCRETE NEWTON METHOD WITH DIRECT SOLUTION OF THE
*         TRUST-REGION SUBPROBLEM.
*  S   PFSED3  COMPRESSED SPARSE STRUCTURE OF THE HESSIAN MATRIX IS
*         COMPUTED FROM THE COORDINATE FORM.
*  S   MXVICP COPYING OF AN INTEGER VECTOR.
*
* EXTERNAL SUBROUTINES :
*  SE  OBJ  COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION.
*         CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS THE NUMBER
*         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND FF IS
*         THE VALUE OF THE OBJECTIVE FUNCTION.
*  SE  DOBJ  COMPUTATION OF THE GRADIENT OF THE OBJECTIVE FUNCTION.
*         CALLING SEQUENCE: CALL DOBJ(NF,X,GF) WHERE NF IS THE NUMBER
*         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND GF(NF)
*         IS THE GRADIENT OF THE OBJECTIVE FUNCTION.
*
      SUBROUTINE PNEDS(NF,MH,X,IX,XL,XU,IH,JH,IPAR,RPAR,F,GMAX,ISPAS,
     & IPRNT,ITERM)
      INTEGER NF,MH,IX(*),IH(*),JH(*),IPAR(7),ISPAS,IPRNT,ITERM
      DOUBLE PRECISION X(*),XL(*),XU(*),RPAR(9),F,GMAX
      INTEGER NB,LGF,LHF,LS,LXO,LGO,LXS,LCOL,LPSL,LPERM,
     & LINVP,LWN11,LWN12,LWN13,LWN14,LJH,IFIL,IER
      INTEGER NRES,NDEC,NIN,NIT,NFV,NFG,NFH
      COMMON /STAT/ NRES,NDEC,NIN,NIT,NFV,NFG,NFH
      INTEGER IA(:)
      DOUBLE PRECISION RA(:)
      ALLOCATABLE IA,RA
      IFIL=IPAR(7)
      IF (IFIL.LE.0) IFIL=1
      IF (ISPAS.LE.1) THEN
      CALL PFSED3(NF,MH,IH,JH,IER)
      IF (IER.NE.0) THEN
      WRITE (6,'(''INPUT ERROR : IER = '',I3)') IER
      STOP
      END IF
      ELSE
      MH=IH(NF+1)-1
      END IF
      ALLOCATE (IA(8*NF+5+(2*IFIL+3)*MH),RA(5*NF+1+(2*IFIL+3)*MH))
      NB=1
*
*     POINTERS FOR AUXILIARY ARRAYS
*
      LGF=1
      LS=LGF+NF
      LXO=LS+NF
      LGO=LXO+NF
      LXS=LGO+NF+1
      LHF=LXS+NF
      LCOL=1
      LPSL=LCOL+NF
      LPERM=LPSL+NF+1
      LINVP=LPERM+NF
      LWN11=LINVP+NF
      LWN12=LWN11+NF+1
      LWN13=LWN12+NF+1
      LWN14=LWN13+NF+1
      LJH=LWN14+NF+1
      CALL MXVICP(IH(NF+1)-1,JH,IA(LJH))
      CALL PNED(NF,NB,(2*IFIL+3)*MH,X,IX,XL,XU,RA(LGF),RA(LHF),IH,
     & IA(LJH),RA(LS),RA(LXO),RA(LGO),RA(LXS),IA(LCOL),IA(LPSL),
     & IA(LPERM),IA(LINVP),IA(LWN11),IA(LWN12),IA(LWN13),IA(LWN14),
     & RPAR(1),RPAR(2),RPAR(3),RPAR(4),RPAR(5),RPAR(6),RPAR(7),GMAX,
     & F,IPAR(1),IPAR(2),IPAR(3),IPAR(4),IPAR(5),IPRNT,ITERM)
      DEALLOCATE (IA,RA)
      RETURN
      END
************************************************************************
* SUBROUTINE PNED               ALL SYSTEMS                   01/09/22
* PURPOSE :
* GENERAL SUBROUTINE FOR LARGE-SCALE BOX CONSTRAINED MINIMIZATION OF
* FUNCTIONS WITH SPARSE HESSIAN MATRICES BASED ON THE DISCRETE NEWTON
* METHOD WITH DIRECT SOLUTION OF THE TRUST-REGION SUBPROBLEM.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  II  NB  CHOICE OF SIMPLE BOUNDS. NB=0-SIMPLE BOUNDS SUPPRESSED.
*         NB>0-SIMPLE BOUNDS ACCEPTED.
*  II  MMAX  MAXIMUM DIMENSION OF THE SPARSE TABLEAU.
*  RU  X(NF)  VECTOR OF VARIABLES.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS. IX(I)=0-VARIABLE
*         X(I) IS UNBOUNDED. IX(I)=1-LOWER BOUND XL(I).LE.X(I).
*         IX(I)=2-UPPER BOUND X(I).LE.XU(I). IX(I)=3-TWO SIDE BOUND
*         XL(I).LE.X(I).LE.XU(I). IX(I)=5-VARIABLE X(I) IS FIXED.
*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
*  RU  GF(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RA  HF(MMAX)  NONZERO ELEMENTS OF THE SPARSE HESSIAN MATRIX
*         TOGETHER WITH AN ADDITIONAL SPACE USED FOR THE NUMERICAL
*         DIFFERENTIATION.
*  II  IH(NF+1)  POINTERS OF DIAGONAL ELEMENTS OF THE MATRIX H.
*  IU  JH(MMAX)  INDICES OF NONZERO ELEMENTS OF THE MATRIX H
*         TOGETHER WITH AN ADDITIONAL SPACE USED FOR THE NUMERICAL
*         DIFFERENTIATION.
*  RU  S(NF)  DIRECTION VECTOR.
*  RU  XO(NF)  VECTORS OF VARIABLES DIFFERENCE.
*  RU  GO(NF)  GRADIENTS DIFFERENCE.
*  RA  XS(NF)  AUXILIARY VECTOR.
*  IA  COL(NF)  AUXILIARY ARRAY.
*  IA  PSL(NF+1)  POINTER VECTOR OF THE COMPACT FORM OF THE TRIANGULAR
*         FACTOR OF THE HESSIAN APPROXIMATION.
*  IA  PERM(NF)  PERMUTATION VECTOR.
*  IA  INVP(NF)  INVERSE PERMUTATION VECTOR.
*  IA  WN11(NF+1) AUXILIARY VECTOR.
*  IA  WN12(NF+1) AUXILIARY VECTOR.
*  IA  WN13(NF+1) AUXILIARY VECTOR.
*  IA  WN14(NF+1) AUXILIARY VECTOR.
*  RI  XMAX  MAXIMUM STEPSIZE.
*  RI  TOLX  TOLERANCE FOR CHANGE OF VARIABLES.
*  RI  TOLF  TOLERANCE FOR CHANGE OF FUNCTION VALUES.
*  RI  TOLB  TOLERANCE FOR THE FUNCTION VALUE.
*  RI  TOLG  TOLERANCE FOR THE GRADIENT NORM.
*  RI  FMIN  ESTIMATION OF THE MINIMUM FUNCTION VALUE.
*  RI  XDEL  TRUST REGION STEPSIZE.
*  RO  GMAX  MAXIMUM PARTIAL DERIVATIVE.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  II  MIT  MAXIMUM NUMBER OF ITERATIONS.
*  II  MFV  MAXIMUM NUMBER OF FUNCTION EVALUATIONS.
*  II  MFG  MAXIMUM NUMBER OF GRADIENT EVALUATIONS.
*  II  IEST  ESTIMATION INDICATOR. IEST=0-MINIMUM IS NOT ESTIMATED.
*         IEST=1-MINIMUM IS ESTIMATED BY THE VALUE FMIN.
*  II  MOS  METHOD FOR COMPUTING THE TRUST REGION STEP. MOS=1-THE
*         DOG-LEG METHOD. MOS=2-THE MORE-SORENSEN METHOD.
*  II  IPRNT  PRINT SPECIFICATION. IPRNT=0-NO PRINT.
*         ABS(IPRNT)=1-PRINT OF FINAL RESULTS.
*         ABS(IPRNT)=2-PRINT OF FINAL RESULTS AND ITERATIONS.
*         IPRNT>0-BASIC FINAL RESULTS. IPRNT<0-EXTENDED FINAL
*         RESULTS.
*  IO  ITERM  VARIABLE THAT INDICATES THE CAUSE OF TERMINATION.
*         ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN
*                   MTESX (USUALLY TWO) SUBSEQUENT ITERATIONS.
*         ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN
*                   MTESF (USUALLY TWO) SUBSEQUENT ITERATIONS.
*         ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB.
*         ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG.
*         ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED,
*                   BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE.
*         ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV.
*         ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED.
*         VALUES ITERM<=-40 DETECT A LACK OF SPACE.
*
* VARIABLES IN COMMON /STAT/ (STATISTICS) :
*  IO  NRES  NUMBER OF RESTARTS.
*  IO  NDEC  NUMBER OF MATRIX DECOMPOSITIONS.
*  IO  NIN  NUMBER OF INNER ITERATIONS.
*  IO  NIT  NUMBER OF ITERATIONS.
*  IO  NFV  NUMBER OF FUNCTION EVALUATIONS.
*  IO  NFG  NUMBER OF GRADIENT EVALUATIONS.
*  IO  NFH  NUMBER OF HESSIAN EVALUATIONS.
*
* SUBPROGRAMS USED :
*  S   PCBS04  ELIMINATION OF BOX CONSTRAINT VIOLATIONS.
*  S   PDSGM1  DIRECTION DETERMINATION USING THE DOUBLE DOG-LEG
*         TRUST-REGION METHOD.
*  S   PDSGM7  DIRECTION DETERMINATION USING THE MORE-SORENSEN
*         TRUST-REGION METHOD.
*  S   PF1HS2  NUMERICAL COMPUTATION OF THE HESSIAN MATRIX USING
*         DIFFERENCES OF GRADIENTS.
*  S   PS0G01  STEPSIZE SELECTION USING TRUST REGION.
*  S   PYADC0  ADDITION OF A BOX CONSTRAINT.
*  S   PYFUT1  TEST ON TERMINATION.
*  S   PYPTSH  DETERMINATION OF GROUPS FOR NUMERICAL DIFFERENTIATION.
*  S   PYRMC0  DELETION OF A BOX CONSTRAINT.
*  S   PYTRCD  COMPUTATION OF PROJECTED DIFFERENCES FOR THE VARIABLE
*         METRIC UPDATE.
*  S   PYTRCG  COMPUTATION OF THE PROJECTED GRADIENT.
*  S   PYTRCS  COMPUTATION OF THE PROJECTED DIRECTION VECTOR.
*  S   PYTSCH  CORRECTION OF THE HESSIAN MATRIX.
*  S   MXSPCC  SPARSE MATRIX REORDERING, SYMBOLIC FACTORIZATION, DATA
*         STRUCTURES TRANSFORMATION. INITIATION OF THE DIRECT SPARSE
*         SOLVER.
*  S   MXSSMI  SPARSE SYMMETRIC MATRIX IS REPLACED BY THE UNIT MATRIX.
*  S   MXUDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVINB  PROJECTION OF A SPARSE SYMMETRIC MATRIX TO SATISFY BOX
*         CONSTRAINTS.
*  S   MXVINE  RESTORATION OF A SPARSE SYMMETRIC MATRIX OBTAINED BY
*         MXVINB
*
* EXTERNAL SUBROUTINES :
*  SE  OBJ  COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION.
*         CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS THE NUMBER
*         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND FF IS
*         THE VALUE OF THE OBJECTIVE FUNCTION.
*  SE  DOBJ  COMPUTATION OF THE GRADIENT OF THE OBJECTIVE FUNCTION.
*         CALLING SEQUENCE: CALL DOBJ(NF,X,GF) WHERE NF IS THE NUMBER
*         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND GF(NF)
*         IS THE GRADIENT OF THE OBJECTIVE FUNCTION.
*
* METHOD :
* DISCRETE NEWTON METHOD WITH TRUST-REGION STRATEGIES BASED ON DIRECT
* MATRIX DECOMPOSITIONS.
*
      SUBROUTINE PNED(NF,NB,MMAX,X,IX,XL,XU,GF,HF,IH,JH,S,XO,GO,XS,
     & COL,PSL,PERM,INVP,WN11,WN12,WN13,WN14,XMAX,TOLX,TOLF,TOLB,
     & TOLG,FMIN,XDEL,GMAX,F,MIT,MFV,MFG,IEST,MOS,IPRNT,ITERM)
      INTEGER NF,NB,MMAX,IX(*),IH(*),JH(*),COL(*),PSL(*),PERM(*),
     & INVP(*),WN11(*),WN12(*),WN13(*),WN14(*),MIT,MFV,MFG,IEST,MOS,
     & IPRNT,ITERM
      DOUBLE PRECISION X(*),XL(*),XU(*),GF(*),HF(*),S(*),XO(*),GO(*),
     & XS(*),TOLX,TOLF,TOLG,TOLB,XDEL,FMIN,XMAX,GMAX,F
      INTEGER IDECF,ITERD,ITERS,ITERH,KD,LD,NTESX,NTESF,MTESX,MTESF,
     & MRED,KIT,IREST,KBF,MES1,MES2,MES3,MAXST,IDIR,ISYS,ITES,KTERS,
     & IRES1,IRES2,NRED,INEW,IOLD,I,M,MH,N
      DOUBLE PRECISION R,RO,FO,FP,P,PO,PP,GNORM,SNORM,RMAX,UMAX,FMAX,
     & DMAX,ETA0,ETA2,ETA9,EPS4,EPS5,EPS8,EPS9,ALF2,BET1,BET2,GAM1,
     & GAM2,DEL1,DEL2,XDELO
      INTEGER NRES,NDEC,NIN,NIT,NFV,NFG,NFH
      COMMON /STAT/ NRES,NDEC,NIN,NIT,NFV,NFG,NFH
      IF (ABS(IPRNT).GT.1) WRITE(6,'(1X,''ENTRY TO PNED :'')')
*
*     INITIATION
*
      KBF=0
      IF (NB.GT.0) KBF=2
      NRES=0
      NDEC=0
      NIN=0
      NIT=0
      NFV=0
      NFG=0
      NFH=0
      IDIR=0
      ISYS=0
      ITES=1
      MTESX=2
      MTESF=2
      ITERM=0
      ITERD=0
      ITERS=2
      ITERH=0
      KTERS=0
      IREST=0
      IRES1=999
      IRES2=0
      IDECF=0
      MRED=10
      MES1=3
      MES2=2
      MES3=1
      ETA0=1.0D-15
      ETA2=1.0D-18
      ETA9=1.0D 120
      EPS4=0.10D 0
      EPS5=0.90D 0
      EPS8=1.00D 0
      EPS9=1.00D-8
      ALF2=1.0D 6
      BET1=0.05D 0
      BET2=0.75D 0
      GAM1=2.0D 0
      GAM2=1.0D 6
      DEL1=0.9D 0
      DEL2=1.1D 0
      RMAX=ETA9
      DMAX=ETA9
      FMAX=1.0D 60
      IF (IEST.LE.0) FMIN=-1.0D 60
      IF (IEST.GT.0) IEST=1
      IF (XMAX.LE.0.0D 0) XMAX=1.0D 16
      XDEL=MIN(XDEL,XMAX)
      IF (TOLX.LE.0.0D 0) TOLX=1.0D-16
      IF (TOLF.LE.0.0D 0) TOLF=1.0D-14
      IF (TOLG.LE.0.0D 0) TOLG=1.0D-6
      IF (TOLB.LE.0.0D 0) TOLB=FMIN+1.0D-16
      IF (MIT.LE.0) MIT=5000
      IF (MFV.LE.0) MFV=5000
      IF (MFG.LE.0) MFG=10000
      IF (MOS.LE.0) MOS=2
      KD= 2
      LD=-1
      KIT=0
      FO=FMIN
*
*     INITIAL OPERATIONS WITH SIMPLE BOUNDS
*
      IF (KBF.GT.0) THEN
      DO 2 I = 1,NF
      IF ((IX(I).EQ.3.OR.IX(I).EQ.4) .AND. XU(I).LE.XL(I)) THEN
      XU(I) = XL(I)
      IX(I) = 5
      ELSE IF (IX(I).EQ.5 .OR. IX(I).EQ.6) THEN
      XL(I) = X(I)
      XU(I) = X(I)
      IX(I) = 5
      END IF
    2 CONTINUE
      END IF
      M=IH(NF+1)-1
      IF (KBF.GT.0) THEN
      CALL PCBS04(NF,X,IX,XL,XU,EPS9,KBF)
      CALL PYADC0(NF,N,X,IX,XL,XU,INEW)
      CALL MXVINE(IH(NF+1)-1,JH)
      END IF
      CALL PYPTSH(NF,MMAX,IH,JH,COL,S,XO,GO,WN11,WN12,GF,ITERM)
      MH=0
      CALL MXSPCC(NF,M,MH,MMAX,HF,IH,JH,PSL,PERM,INVP,WN11,WN12,WN13,
     & WN14,ITERM)
      IF (ITERM.NE.0) GO TO 11180
      CALL OBJ(NF,X,F)
      NFV=NFV+1
      CALL DOBJ(NF,X,GF)
      NFG=NFG+1
12140 CONTINUE
      CALL PF1HS2(NF,MH,MMAX,X,IX,S,HF,IH,JH,GO,GF,COL,WN11,WN12,XS,
     & F,ETA0,KBF,ITERM,ISYS)
      IF (ISYS.EQ.0) GO TO 12160
      CALL DOBJ(NF,X,GO)
      NFG=NFG+1
      GO TO 12140
12160 CONTINUE
      IDECF=0
*
*     START OF THE ITERATION WITH TESTS FOR TERMINATION.
*
11120 CONTINUE
      CALL PYTRCG(NF,NF,IX,GF,UMAX,GMAX,KBF,IOLD)
      IF (ABS(IPRNT).GT.1)
     & WRITE (6,'(1X,''NIT='',I5,2X,''NFV='',I5,2X,''NFG='',I5,2X,
     & ''F='', G16.9,2X,''G='',E10.3)') NIT,NFV,NFG,F,GMAX
      CALL PYFUT1(NF,F,FO,UMAX,GMAX,DMAX,TOLX,TOLF,TOLB,TOLG,KD,
     & NIT,KIT,MIT,NFV,MFV,NFG,MFG,NTESX,MTESX,NTESF,MTESF,ITES,
     & IRES1,IRES2,IREST,ITERS,ITERM)
      IF (ITERM.NE.0) GO TO 11180
      IF (KBF.GT.0) CALL PYRMC0(NF,N,IX,GF,EPS8,UMAX,GMAX,RMAX,IOLD,
     & IREST)
11130 CONTINUE
      IF (IREST.GT.0) THEN
      CALL MXSSMI(NF,HF,IH)
      IDECF=0
      LD=MIN(LD,1)
      IF (KIT.LT.NIT) THEN
        NRES=NRES+1
        KIT = NIT
      ELSE
        ITERM=-10
        IF (ITERS.LT.0) ITERM=ITERS-5
        IF (GMAX.LE.1.0D 2*TOLG) ITERM=-ITERM
      END IF
      END IF
      IF (ITERM.NE.0) GO TO 11180
      IF (KBF.GT.0) THEN
      CALL MXVINB(M,IX,JH)
      CALL PYTSCH(NF,IX,HF,IH,JH,KBF)
      END IF
*
*     DIRECTION DETERMINATION
*
      IF (MOS.LE.1) THEN
      CALL PDSGM1(NF,MMAX,MH,IX,GF,HF,IH,JH,S,XO,GO,XS,PSL,PERM,WN11,
     & WN12,XMAX,XDEL,GNORM,SNORM,FMIN,F,P,PP,ETA2,ALF2,KD,KBF,IEST,
     & IDECF,NDEC,ITERD,ITERM)
      ELSE
      CALL PDSGM7(NF,MMAX,MH,IX,GF,HF,IH,JH,S,XO,GO,PSL,PERM,WN11,WN12,
     & XMAX,XDEL,XDELO,GNORM,SNORM,FMIN,F,P,PP,ETA2,DEL1,DEL2,KD,KBF,
     & IEST,IDIR,IDECF,NDEC,ITERD,ITERM)
      END IF
*
*     TEST ON LOCALLY CONSTRAINED STEP AND PREPARATION OF STEPSIZE
*     SELECTION
*
      IF (ITERD.LT.0) THEN
        ITERM=ITERD
      ELSE
      IF (SNORM.LE.0.0D 0) THEN
        IREST=MAX(IREST,1)
      ELSE
        IREST=0
      END IF
      IF (IREST.EQ.0) THEN
        RMAX=XMAX/SNORM
      END IF
      END IF
      IF (ITERM.NE.0) GO TO 11180
      IF (IREST.NE.0) GO TO 11130
      IF (NIT.EQ.1) KIT=NIT
      CALL PYTRCS(NF,X,IX,XO,XL,XU,GF,GO,S,RO,FP,FO,F,PO,P,RMAX,ETA9,
     & KBF)
      IF (RMAX.EQ.0.0D 0) GO TO 11175
11150 CONTINUE
      CALL PS0G01(R,F,FO,PO,PP,XDEL,XDELO,XMAX,RMAX,SNORM,
     & BET1,BET2,GAM1,GAM2,EPS4,EPS5,KD,LD,IDIR,ITERS,ITERD,
     & MAXST,NRED,MRED,KTERS,MES1,MES2,MES3,ISYS)
      IF (ISYS.EQ.0) GO TO 11154
      CALL MXUDIR(NF,R,S,XO,X,IX,KBF)
      CALL PCBS04(NF,X,IX,XL,XU,EPS9,KBF)
      CALL OBJ(NF,X,F)
      NFV=NFV+1
      GO TO 11150
11154 CONTINUE
      KD=2
      IF (ITERS.LE.0) THEN
      R=0.0D 0
      F=FO
      P=PO
      CALL MXVCOP(NF,XO,X)
      IF (ITERS.LT.0) THEN
      ITERM=-6
      IF (GMAX.LE.1.0D 2*TOLG) ITERM=-ITERM
      GO TO 11180
      END IF
      IF (IDIR.EQ.0) IREST=MAX(IREST,1)
      LD=KD
      GO TO 11130
      END IF
      IF (KD.GT.LD) THEN
      CALL DOBJ(NF,X,GF)
      NFG=NFG+1
13140 CONTINUE
      CALL PF1HS2(NF,MH,MMAX,X,IX,S,HF,IH,JH,GO,GF,COL,WN11,WN12,XS,
     & F,ETA0,KBF,ITERM,ISYS)
      IF (ISYS.EQ.0) GO TO 13160
      CALL DOBJ(NF,X,GO)
      NFG=NFG+1
      GO TO 13140
13160 CONTINUE
      IDECF=0
      END IF
      CALL PYTRCD(NF,X,IX,XO,GF,GO,R,F,FO,P,PO,DMAX,KBF,KD,LD,ITERS)
11175 CONTINUE
      IF (KBF.GT.0) CALL PYADC0(NF,N,X,IX,XL,XU,INEW)
      GO TO 11120
11180 CONTINUE
      IF (IPRNT.GT.1.OR.IPRNT.LT.0)
     & WRITE(6,'(1X,''EXIT FROM PNED :'')')
      IF (IPRNT.NE.0)
     & WRITE (6,'(1X,''NIT='',I5,2X,''NFV='',I5,2X,''NFG='',I5,2X,
     & ''F='', G16.9,2X,''G='',E10.3,2X,''ITERM='',I3)') NIT,NFV,NFG,
     & F,GMAX,ITERM
      IF (IPRNT.LT.0)
     & WRITE (6,'(1X,''X='',5(G14.7,1X):/(3X,5(G14.7,1X)))')
     & (X(I),I=1,NF)
      RETURN
      END