************************************************************************
* SUBROUTINE PGACU              ALL SYSTEMS                   97/01/22
* PURPOSE :
* EASY TO USE SUBROUTINE FOR LARGE-SCALE UNCONSTRAINED LEAST-SQUARES
* PROBLEMS WITH SPARSE JACOBIAN MATRICES.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  II  NA  NUMBER OF APPROXIMATED FUNCTIONS.
*  IU  MA  NUMBER OF NONZERO ELEMENTS IN THE JACOBIAN MATRIX.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  RO  AF(NA)  VECTOR CONTAINING VALUES OF THE APPROXIMATED
*         FUNCTIONS.
*  RI  IAG(NA+1)  POSITION OF THE FIRST ROWS ELEMENTS IN THE JACOBIAN
*         MATRIX.
*  RI  JAG(MA)  COLUMN INDICES OF ELEMENTS IN THE JACOBIAN MATRIX.
*  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)  TYPE OF THE SECOND ORDER CORRECTION. MET=1-THE MARWIL
*         SPARSE VARIABLE METRIC UPDATE. MET=2-THE NEWTON CORRECTION
*         BY USING GRADIENT DIFFERENCES. MET=3-THE PARTITIONED VARIABLE
*         METRIC UPDATE.
*      IPAR(5)  METHOD FOR COMPUTING A TRUST REGION STEP. IPAR(5)=1-THE
*         STEIHAUG-TOINT METHOD. IPAR(5)=2-THE SHIFTED STEIHAUG-TOINT
*         METHOD WITH FIVE LANCZOS STEPS. IPAR(5)>2-THE SHIFTED
*         STEIHAUG-TOINT METHOD WITH IPAR(5) LANCZOS STEPS.
*      IPAR(6)  TYPE OF PRECONDITIONING. IPAR(6)=1-PRECONDITIONING IS
*         NOT USED. IPAR(6)=2-PRECONDITIONING BY THE INCOMPLETE
*         GILL-MURRAY DECOMPOSITION. IPAR(6)=3-PRECONDITIONING BY THE
*         INCOMPLETE GILL-MURRAY DECOMPOSITION WITH A PRELIMINARY
*         SOLUTION OF THE PRECONDITIONED SYSTEM WHICH IS USED IF IT
*         SATISFIES THE TERMINATION CRITERION.
*      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)  TOLERANCE FOR THE SWITCH BETWEEN THE GAUSS-NEWTON AND
*      RPAR(9)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PNED.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RO  GMAX  MAXIMUM PARTIAL DERIVATIVE.
*  II  IDER  DEGREE OF ANALYTICALLY COMPUTED DERIVATIVES (0 OR 1).
*  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   PGAC  HYBRID METHOD FOR SPARSE LEAST SQUARES PROBLEMS.
*  S   PASED3  COMPRESSED SPARSE STRUCTURE OF THE JACOBIAN MATRIX IS
*         COMPUTED FROM THE COORDINATE FORM.
*  S   PFSET2  NUMBER OF NONZERO ELEMENTS IN THE PARTITIONED HESSIAN
*         MATRIX.
*
* EXTERNAL SUBROUTINES :
*  SE  FUN  COMPUTATION OF THE VALUE OF THE APPROXIMATED FUNCTION.
*         CALLING SEQUENCE: CALL FUN(NF,KA,X,FA) WHERE NF IS A NUMBER
*         OF VARIABLES, KA IS THE INDEX OF THE APPROXIMATED FUNCTION,
*         X(NF) IS A VECTOR OF VARIABLES AND FA IS THE VALUE OF THE
*         APPROXIMATED FUNCTION.
*  SE  DFUN  COMPUTATION OF THE GRADIENT OF THE APPROXIMATED FUNCTION.
*         CALLING SEQUENCE: CALL DFUN(NF,KA,X,GA) WHERE NF IS A NUMBER
*         OF VARIABLES, KA IS THE INDEX OF THE APPROXIMATED FUNCTION,
*         X(NF) IS A VECTOR OF VARIABLES AND GA(NF) IS THE GRADIENT OF
*         THE APPROXIMATED FUNCTION.
*
      SUBROUTINE PGACU(NF,NA,MA,X,AF,IAG,JAG,IPAR,RPAR,F,GMAX,IDER,
     & ISPAS,IPRNT,ITERM)
      INTEGER NF,NA,MA,IAG(*),JAG(*),IPAR(7),IDER,ISPAS,IPRNT,ITERM
      DOUBLE PRECISION X(*),AF(*),RPAR(9),F,GMAX
      INTEGER NB,LGA,LAG,LG,LHA,LAH,LH,LS,LXO,LGO,LAGO,LXS,LGS,LIH,LJH,
     & LIW,ML,MH,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 PASED3(NF,NA,MA,IAG,JAG,IER)
      IF (IER.NE.0) THEN
      WRITE (6,'(''INPUT ERROR : IER = '',I3)') IER
      STOP
      END IF
      ELSE
      MA=IAG(NA+1)-1
      END IF
      CALL PFSET2(NA,MH,ML,IAG)
      ALLOCATE (IA(2*NF+1+(IFIL+2)*MH))
      IF (IPAR(4).LE.2) THEN
      ALLOCATE (RA(7*NF+ML+(IFIL+2)*MH))
      ELSE
      ALLOCATE (RA(7*NF+2*MA+ML+(IFIL+3)*MH))
      END IF
      NB=0
*
*     POINTERS FOR AUXILIARY ARRAYS
*
      LGA=1
      LG=LGA+NF
      LS=LG+NF
      LXO=LS+NF
      LGO=LXO+NF
      LXS=LGO+NF
      LGS=LXS+NF
      IF (IPAR(4).LE.2) THEN
      LAG=LGS
      LAGO=LGS
      LAH=LGS
      LHA=LGS+NF
      ELSE
      LAG=LGS+NF
      LAGO=LAG+MA
      LAH=LAGO+MA
      LHA=LAH+MH
      END IF
      LH=LHA+ML
      LIH=1
      LIW=LIH+NF+1
      LJH=LIW+NF+1
      CALL PGAC(NF,NA,NB,(IFIL+2)*MH,X,IA,RA,RA,AF,RA(LGA),RA(LAG),
     & RA(LG),RA(LHA),RA(LAH),RA(LH),IA(LIH),IA(LJH),IAG,JAG,RA(LS),
     & RA(LXO),RA(LGO),RA(LAGO),RA(LXS),RA(LGS),IA(LIW),RPAR(1),
     & RPAR(2),RPAR(3),RPAR(4),RPAR(5),RPAR(6),RPAR(7),RPAR(8),GMAX,
     & F,IPAR(1),IPAR(2),IPAR(3),IPAR(4),IPAR(5),IPAR(6),IDER,IPRNT,
     & ITERM)
      DEALLOCATE (IA,RA)
      RETURN
      END
************************************************************************
* SUBROUTINE PGACS              ALL SYSTEMS                   97/01/22
* PURPOSE :
* EASY TO USE SUBROUTINE FOR LARGE-SCALE BOX CONSTRAINED LEAST-SQUARES
* PROBLEMS WITH SPARSE JACOBIAN MATRICES.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  II  NA  NUMBER OF LINEAR APPROXIMATED FUNCTIONS.
*  IU  MA  NUMBER OF NONZERO ELEMENTS IN THE JACOBIAN 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.
*  RO  AF(NA)  VECTOR CONTAINING VALUES OF THE APPROXIMATED
*         FUNCTIONS.
*  RI  IAG(NA+1)  POSITION OF THE FIRST ROWS ELEMENTS IN THE JACOBIAN
*         MATRIX.
*  RI  JAG(MA)  COLUMN INDICES OF ELEMENTS IN THE JACOBIAN MATRIX.
*  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)  TYPE OF THE SECOND ORDER CORRECTION. MET=1-THE MARWIL
*         SPARSE VARIABLE METRIC UPDATE. MET=2-THE NEWTON CORRECTION
*         BY USING GRADIENT DIFFERENCES. MET=3-THE PARTITIONED VARIABLE
*         METRIC UPDATE.
*      IPAR(5)  METHOD FOR COMPUTING A TRUST REGION STEP. IPAR(5)=1-THE
*         STEIHAUG-TOINT METHOD. IPAR(5)=2-THE SHIFTED STEIHAUG-TOINT
*         METHOD WITH FIVE LANCZOS STEPS. IPAR(5)>2-THE SHIFTED
*         STEIHAUG-TOINT METHOD WITH IPAR(5) LANCZOS STEPS.
*      IPAR(6)  TYPE OF PRECONDITIONING. IPAR(6)=1-PRECONDITIONING IS
*         NOT USED. IPAR(6)=2-PRECONDITIONING BY THE INCOMPLETE
*         GILL-MURRAY DECOMPOSITION. IPAR(6)=3-PRECONDITIONING BY THE
*         INCOMPLETE GILL-MURRAY DECOMPOSITION WITH A PRELIMINARY
*         SOLUTION OF THE PRECONDITIONED SYSTEM WHICH IS USED IF IT
*         SATISFIES THE TERMINATION CRITERION.
*      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)  TOLERANCE FOR THE SWITCH BETWEEN THE GAUSS-NEWTON AND
*      RPAR(9)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PNED.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RO  GMAX  MAXIMUM PARTIAL DERIVATIVE.
*  II  IDER  DEGREE OF ANALYTICALLY COMPUTED DERIVATIVES (0 OR 1).
*  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   PGAC  HYBRID METHOD FOR SPARSE LEAST SQUARES PROBLEMS.
*  S   PASED3  COMPRESSED SPARSE STRUCTURE OF THE JACOBIAN MATRIX IS
*         COMPUTED FROM THE COORDINATE FORM.
*  S   PFSET2  NUMBER OF NONZERO ELEMENTS IN THE PARTITIONED HESSIAN
*         MATRIX.
*
* EXTERNAL SUBROUTINES :
*  SE  FUN  COMPUTATION OF THE VALUE OF THE APPROXIMATED FUNCTION.
*         CALLING SEQUENCE: CALL FUN(NF,KA,X,FA) WHERE NF IS A NUMBER
*         OF VARIABLES, KA IS THE INDEX OF THE APPROXIMATED FUNCTION,
*         X(NF) IS A VECTOR OF VARIABLES AND FA IS THE VALUE OF THE
*         APPROXIMATED FUNCTION.
*  SE  DFUN  COMPUTATION OF THE GRADIENT OF THE APPROXIMATED FUNCTION.
*         CALLING SEQUENCE: CALL DFUN(NF,KA,X,GA) WHERE NF IS A NUMBER
*         OF VARIABLES, KA IS THE INDEX OF THE APPROXIMATED FUNCTION,
*         X(NF) IS A VECTOR OF VARIABLES AND GA(NF) IS THE GRADIENT OF
*         THE APPROXIMATED FUNCTION.
*
      SUBROUTINE PGACS(NF,NA,MA,X,IX,XL,XU,AF,IAG,JAG,IPAR,RPAR,F,GMAX,
     & IDER,ISPAS,IPRNT,ITERM)
      INTEGER NF,NA,MA,IX(*),IAG(*),JAG(*),IPAR(7),IDER,ISPAS,IPRNT,
     & ITERM
      DOUBLE PRECISION X(*),XL(*),XU(*),AF(*),RPAR(9),F,GMAX
      INTEGER NB,LGA,LG,LAG,LHA,LAH,LH,LS,LXO,LGO,LAGO,LXS,LGS,LIH,LJH,
     & LIW,ML,MH,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 PASED3(NF,NA,MA,IAG,JAG,IER)
      IF (IER.NE.0) THEN
      WRITE (6,'(''INPUT ERROR : IER = '',I3)') IER
      STOP
      END IF
      ELSE
      MA=IAG(NA+1)-1
      END IF
      CALL PFSET2(NA,MH,ML,IAG)
      ALLOCATE (IA(2*NF+1+(IFIL+2)*MH))
      IF (IPAR(4).LE.2) THEN
      ALLOCATE (RA(7*NF+ML+(IFIL+2)*MH))
      ELSE
      ALLOCATE (RA(7*NF+2*MA+ML+(IFIL+3)*MH))
      END IF
      NB=1
*
*     POINTERS FOR AUXILIARY ARRAYS
*
      LGA=1
      LG=LGA+NF
      LS=LG+NF
      LXO=LS+NF
      LGO=LXO+NF
      LXS=LGO+NF
      LGS=LXS+NF
      IF (IPAR(4).LE.2) THEN
      LAG=LGS
      LAGO=LGS
      LAH=LGS
      LHA=LGS+NF
      ELSE
      LAG=LGS+NF
      LAGO=LAG+MA
      LAH=LAGO+MA
      LHA=LAH+MH
      END IF
      LH=LHA+ML
      LIH=1
      LIW=LIH+NF+1
      LJH=LIW+NF+1
      CALL PGAC(NF,NA,NB,(IFIL+2)*MH,X,IX,XL,XU,AF,RA(LGA),RA(LAG),
     & RA(LG),RA(LHA),RA(LAH),RA(LH),IA(LIH),IA(LJH),IAG,JAG,RA(LS),
     & RA(LXO),RA(LGO),RA(LAGO),RA(LXS),RA(LGS),IA(LIW),RPAR(1),
     & RPAR(2),RPAR(3),RPAR(4),RPAR(5),RPAR(6),RPAR(7),RPAR(8),GMAX,
     & F,IPAR(1),IPAR(2),IPAR(3),IPAR(4),IPAR(5),IPAR(6),IDER,IPRNT,
     & ITERM)
      DEALLOCATE (IA,RA)
      RETURN
      END
************************************************************************
* SUBROUTINE PGAC               ALL SYSTEMS                   01/09/22
* PURPOSE :
* GENERAL SUBROUTINE FOR LARGE-SCALE BOX CONSTRAINED LEAST-SQUARES
* PROBLEMS WITH SPARSE JACOBIAN MATRICES.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  II  NA  NUMBER OF APPROXIMATING FUNCTIONS.
*  II  NB  CHOICE OF SIMPLE BOUNDS. NB=0-SIMPLE BOUNDS SUPPRESSED.
*         NB>0-SIMPLE BOUNDS ACCEPTED.
*  II  MMAX  MAXIMUM DIMENSION OF THE SPARSE TABLEAU.
*  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.
*  RO  AF(NA)  VECTOR CONTAINING VALUES OF THE APPROXIMATED
*         FUNCTIONS.
*  RA  GA(NF)  GRADIENT OF THE SELECTED APPROXIMATING FUNCTION.
*  RA  AG(MA)  JACOBIAN MATRIX OF APPROXIMATING FUNCTIONS.
*  RU  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RA  HA(ML)  HESSIAN MATRIX OF THE SELECTED APPROXIMATING FUNCTION.
*  RA  AH(MH)  ELEMENTS OF THE PARTITIONED HESSIAN MATRIX.
*  RA  H(MMAX)  NONZERO ELEMENTS OF THE SPARSE HESSIAN MATRIX 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.
*  RI  IAG(NA+1)  POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD AG.
*  RI  JAG(MA)  COLUMN INDICES OF ELEMENTS IN THE FIELD AG.
*  RU  S(NF)  DIRECTION VECTOR.
*  RU  XO(NF)  VECTORS OF VARIABLES DIFFERENCE.
*  RU  GO(NF)  GRADIENTS DIFFERENCE.
*  RA  AGO(MA)  OLD JACOBIAN MATRIX OF THE PARTITIONED FUNCTION,
*  RA  XS(NF)  AUXILIARY VECTOR.
*  RA  GS(NF)  AUXILIARY VECTOR.
*  IA  IW(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.
*  RI  ETA   TOLERANCE FOR THE SWITCH BETWEEN THE GAUSS-NEWTON AND
*         THE MARWIL METHODS.
*  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  MEC  TYPE OF THE SECOND ORDER CORRECTION. MEC=1-THE MARWIL
*         SPARSE VARIABLE METRIC UPDATE. MEC=2-THE NEWTON CORRECTION
*         BY USING GRADIENT DIFFERENCES. MEC=3-THE PARTITIONED
*         VARIABLE METRIC UPDATE.
*  II  MOS1  METHOD FOR COMPUTING A TRUST REGION STEP. MOS1=1-THE
*         STEIHAUG-TOINT METHOD. MOS1=2-THE SHIFTED STEIHAUG-TOINT
*         METHOD WITH FIVE LANCZOS STEPS. MOS1>2-THE SHIFTED
*         STEIHAUG-TOINT METHOD WITH MOS1 LANCZOS STEPS.
*  II  MOS2  TYPE OF PRECONDITIONING. MOS2=1-PRECONDITIONING IS NOT
*         USED. MOS2=2-PRECONDITIONING BY THE INCOMPLETE GILL-MURRAY
*         DECOMPOSITION. MOS2=3-PRECONDITIONING BY THE INCOMPLETE
*         GILL-MURRAY DECOMPOSITION WITH A PRELIMINARY SOLUTION OF
*         THE PRECONDITIONED SYSTEM WHICH IS USED IF IT SATISFIES
*         THE TERMINATION CRITERION.
*  II  IDER  DEGREE OF ANALYTICALLY COMPUTED DERIVATIVES (0 OR 1).
*  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   PA2SQ4  COMPUTATION OF THE VALUE, THE GRADIENT, AND THE SPARSE
*         NORMAL EQUATION MATRIX OF THE OBJECTIVE FUNCTION WHICH IS
*         DEFINED AS A SUM OF SQUARES OF THE APPROXIMATED FUNCTIONS.
*  S   PA2SQ8  COMPUTATION OF THE VALUE, THE GRADIENT, AND THE SPARSE
*         NORMAL EQUATION MATRIX OF THE OBJECTIVE FUNCTION WHICH IS
*         DEFINED AS A SUM OF SQUARES OF THE APPROXIMATED FUNCTIONS.
*  S   PCBS04  ELIMINATION OF BOX CONSTRAINT VIOLATIONS.
*  S   PDSGM4  DIRECTION DETERMINATION USING THE STEIHAUG-TOINT AND
*         SHIFTED STEIHAUG-TOINT TRUST-REGION METHOD.
*  S   PFSET3  PREPARATION OF THE SPARSE NORMAL EQUATION MATRIX
*         STRUCTURE.
*  S   PS0G01  STEPSIZE SELECTION USING TRUST REGION.
*  S   PUBBM2  PARTITIONED VARIABLE METRIC UPDATES.
*  S   PUSSD5  COMPUTATION OF THE GAUSS-NEWTON MATRIX.
*  S   PUSMM1  SPARSE VARIABLE METRIC (MARWIL) UPDATE.
*  S   PYADC0  ADDITION OF A BOX CONSTRAINT.
*  S   PYFUT1  TEST ON TERMINATION.
*  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   MXBSMI  INITIATION OF THE PARTITIONED HESSIAN MATRIX.
*  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   MXVDIF  DIFFERENCE OF TWO VECTORS.
*  S   MXVINB  PROJECTION OF A SPARSE SYMMETRIC MATRIX TO SATISFY BOX
*         CONSTRAINTS.
*  S   MXVINE  RESTORATION OF A SPARSE SYMMETRIC MATRIX OBTAINED BY
*         MXVINB
*  S   MXVINS  INITIATION OF THE INTEGER VECTOR.
*  S   MXVSAV  RESTORATION OF THE GRADIENT VECTOR.
*
* EXTERNAL SUBROUTINES :
*  SE  FUN  COMPUTATION OF THE VALUE OF THE APPROXIMATED FUNCTION.
*         CALLING SEQUENCE: CALL FUN(NF,KA,X,FA) WHERE NF IS A NUMBER
*         OF VARIABLES, KA IS THE INDEX OF THE APPROXIMATED FUNCTION,
*         X(NF) IS A VECTOR OF VARIABLES AND FA IS THE VALUE OF THE
*         APPROXIMATED FUNCTION.
*  SE  DFUN  COMPUTATION OF THE GRADIENT OF THE APPROXIMATED FUNCTION.
*         CALLING SEQUENCE: CALL DFUN(NF,KA,X,GA) WHERE NF IS A NUMBER
*         OF VARIABLES, KA IS THE INDEX OF THE APPROXIMATED FUNCTION,
*         X(NF) IS A VECTOR OF VARIABLES AND GA(NF) IS THE GRADIENT OF
*         THE APPROXIMATED FUNCTION.
*
* METHOD :
* HYBRID GAUSS-NEWTON METHOD WITH VARIOUS SECOND ORDER CORRECTIONS AND
* TRUST-REGION STRATEGIES BASED ON CONJUGATE GRADIENT ITERATIONS.
*
      SUBROUTINE PGAC(NF,NA,NB,MMAX,X,IX,XL,XU,AF,GA,AG,G,HA,AH,H,IH,
     & JH,IAG,JAG,S,XO,GO,AGO,XS,GS,IW,XMAX,TOLX,TOLF,TOLB,TOLG,FMIN,
     & XDEL,ETA,GMAX,F,MIT,MFV,MFG,MEC,MOS1,MOS2,IDER,IPRNT,ITERM)
      INTEGER NF,NA,NB,MMAX,IX(*),IH(*),JH(*),IAG(*),JAG(*),IW(*),MIT,
     & MFV,MFG,MEC,MOS1,MOS2,IDER,IPRNT,ITERM
      DOUBLE PRECISION X(*),XL(*),XU(*),AF(*),GA(*),AG(*),G(*),HA(*),
     & AH(*),H(*),S(*),XO(*),GO(*),AGO(*),XS(*),GS(*),XMAX,TOLX,TOLF,
     & TOLG,TOLB,XDEL,FMIN,ETA,GMAX,F
      INTEGER IDECF,ITERD,ITERS,ITERH,KD,LD,NTESX,NTESF,MTESX,MTESF,
     & MRED,KIT,IREST,KBF,MET,MET1,MET3,MES1,MES2,MES3,MOS3,MAXST,IDIR,
     & ISYS,IEST,ITES,KTERS,IRES1,IRES2,NRED,IPOM1,INEW,IOLD,I,M,N,ISNA
      DOUBLE PRECISION R,RO,FO,FP,P,PO,PP,GNORM,GNORMO,SNORM,RMAX,FMAX,
     & DMAX,UMAX,ETA0,ETA2,ETA9,EPS4,EPS5,EPS8,EPS9,BET1,BET2,GAM1,
     & GAM2,DEL1,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 PGAC :'')')
*
*     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=-1
      IRES1=999
      IRES2=0
      IDECF=0
      MRED=10
      MES1=3
      MES2=2
      MES3=1
      MOS3=0
      ETA0=1.0D-15
      ETA2=1.0D-8
      ETA9=1.0D 120
      EPS4=0.10D 0
      EPS5=0.90D 0
      EPS8=1.00D 0
      EPS9=1.00D-8
      BET1=0.05D 0
      BET2=0.75D 0
      GAM1=2.0D 0
      GAM2=1.0D 6
      DEL1=0.95D 0
      RMAX=ETA9
      DMAX=ETA9
      FMAX=1.0D 60
      FMIN=MAX(FMIN,0.0D 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 (ETA .LE.0.0D 0) ETA =1.5D-4
      IDER=MAX(IDER,0)
      IF (MIT.LE.0) MIT=5000
      IF (MFV.LE.0) MFV=5000
      IF (MFG.LE.0) MFG=10000
      IF (MEC.LE.0) MEC=2
      IF (MOS1.LE.0) MOS1=2
      IF (MOS1.EQ.2) MOS1=5
      IF (MOS2.LE.0) MOS2=2
      IF (MEC.EQ.1) THEN
      MET1=1
      ELSE IF (MEC.EQ.2) THEN
      MET=2
      ELSE
      MET=4
      MET1=1
      MET3=3
      END IF
      ISNA=MEC-1
      KD= 2
      LD=-1
      KIT=0
      FO=FMIN
      IPOM1=0
*
*     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
      CALL MXVINP(NF+1,IH)
      CALL MXVINP(NF,JH)
      CALL PFSET3(NF,NA,M,MMAX,IH,JH,IAG,JAG,ITERM)
      IF (ITERM.NE.0) GO TO 11180
      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
      IF (ITERM.NE.0) GO TO 11180
      IF (MEC.NE.2) THEN
      CALL PA2SQ4(NF,NA,X,GA,AG,G,H,IH,JH,IAG,JAG,AF,F,ETA0,KD,LD,ISNA,
     & NFV,NFG,IDER,IDECF)
      ELSE
      CALL PA2SQ8(NF,NA,X,IX,GA,G,GO,XS,HA,H,IH,JH,IAG,JAG,AF,F,ETA0,
     & KBF,KD,LD,NFV,NFG,IPOM1,IDER,IDECF)
      END IF
*
*     START OF THE ITERATION WITH TESTS FOR TERMINATION.
*
11120 CONTINUE
      CALL PYTRCG(NF,NF,IX,G,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,G,EPS8,UMAX,GMAX,RMAX,IOLD,
     & IREST)
11130 CONTINUE
      IF (MEC.GE.3) THEN
      IF (IREST.LT.0.OR.NIT.LE.1) CALL MXBSMI(NA,AH,IAG)
      END IF
      IF (IREST.GT.0) THEN
      CALL MXSSMI(NF,H,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
      ELSE IF (MEC.GE.3) THEN
      IF (IDIR.EQ.0.AND.FO-F.LE.ETA*FO) THEN
      CALL PUSSD5(NA,AF,AH,IAG,JAG,H,IH,JH)
      IDECF=0
      LD=MIN(LD,1)
      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,H,IH,JH,KBF)
      END IF
*
*     DIRECTION DETERMINATION
*
      CALL PDSGM4(NF,MMAX,IX,G,H,IH,JH,S,XO,GO,XS,GS,IW,XMAX,
     & XDEL,GNORM,GNORMO,SNORM,FMIN,F,P,PP,ETA0,ETA2,DEL1,KD,KBF,
     & MOS1,MOS2,MOS3,IEST,IDECF,NDEC,NIT,NIN,ITERD,ITERM)
*
*     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,G,GO,S,RO,FP,FO,F,PO,P,RMAX,ETA9,
     & KBF)
      IF (MEC.GE.3) CALL MXVCOP(IAG(NA+1)-1,AG,AGO)
      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)
      IF (MEC.NE.2) THEN
      CALL PA2SQ4(NF,NA,X,GA,AG,G,H,IH,JH,IAG,JAG,AF,F,ETA0,KD,LD,ISNA,
     & NFV,NFG,IDER,IDECF)
      ELSE
      CALL PA2SQ8(NF,NA,X,IX,GA,G,GO,XS,HA,H,IH,JH,IAG,JAG,AF,F,ETA0,
     & KBF,KD,LD,NFV,NFG,IPOM1,IDER,IDECF)
      END IF
      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)
C      IF (MEC.NE.2) THEN
C      ELSE
C      IF (ITERS.LT.0) THEN
C      ITERM=-6
C      IF (GMAX.LE.1.0D 2*TOLG) ITERM=-ITERM
C      GO TO 11180
C      END IF
C      END IF
      IF (IDIR.EQ.0) IREST=MAX(IREST,1)
      LD=KD
      GO TO 11130
      END IF
      IF (MEC.EQ.1) THEN
      IF (FO-F.LE.ETA*FO) KD=1
      ELSE IF (MEC.EQ.2) THEN
      IF (MET.GE.2.AND.FO-F.GE.ETA*FO) IPOM1=1
      END IF
      IF (KD.GT.LD) THEN
      IF (MEC.NE.2) THEN
      CALL PA2SQ4(NF,NA,X,GA,AG,G,H,IH,JH,IAG,JAG,AF,F,ETA0,KD,LD,ISNA,
     & NFV,NFG,IDER,IDECF)
      ELSE
      CALL PA2SQ8(NF,NA,X,IX,GA,G,GO,XS,HA,H,IH,JH,IAG,JAG,AF,F,ETA0,
     & KBF,KD,LD,NFV,NFG,IPOM1,IDER,IDECF)
      END IF
      END IF
      IPOM1=0
      ITERD=0
      CALL PYTRCD(NF,X,IX,XO,G,GO,R,F,FO,P,PO,DMAX,KBF,KD,LD,ITERS)
      IF (MEC.EQ.1) THEN
      IF (FO-F.LE.ETA*FO) THEN
      IDECF=0
      CALL PUSMM1(NF,H,IH,JH,G,XS,S,XO,GO,IX,R,PO,NIT,KIT,MET1,ITERD,
     & ITERH,KBF)
      END IF
      ELSE IF (MEC.EQ.3) THEN
      IF (ITERS.GT.0) THEN
      CALL MXVDIF(IAG(NA+1)-1,AG,AGO,AGO)
      ELSE
      CALL MXVSAV(IAG(NA+1)-1,AG,AGO)
      END IF
      CALL PUBBM2(NA,AH,IAG,JAG,S,XO,AGO,ETA0,ETA9,NIT,KIT,ITERH,MET,
     & MET1,MET3)
      END IF
11175 CONTINUE
      IF (MEC.NE.2) THEN
      IF (IDIR.EQ.0) THEN
      IF (ITERH.NE.0) THEN
      IREST=MAX(IREST,1)
      END IF
      IF (KBF.GT.0) CALL PYADC0(NF,N,X,IX,XL,XU,INEW)
      GO TO 11120
      ELSE
      GO TO 11130
      END IF
      ELSE
      IF (KBF.GT.0) CALL PYADC0(NF,N,X,IX,XL,XU,INEW)
      GO TO 11120
      END IF
11180 CONTINUE
      IF (IPRNT.GT.1.OR.IPRNT.LT.0)
     & WRITE(6,'(1X,''EXIT FROM PGAC :'')')
      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