* SUBROUTINE PYTRND             ALL SYSTEMS                   91/12/01
C PORTABILITY : ALL SYSTEMS
C 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DUAL RANGE SPACE QUADRATIC PROGRAMMING METHOD FOR MINIMAX
* APPROXIMATION.
*
* PARAMETERS :
*  II  NF DECLARED NUMBER OF VARIABLES.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  II  NC NUMBER OF CONSTRAINTS.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  RI  XN(NF)  VECTOR OF SCALING FACTORS.
*  RO  XO(NF)  SAVED VECTOR OF VARIABLES.
*  II  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RO  CZ(NF)  VECTOR OF LAGRANGE MULTIPLIERS.
*  RO  CZS(NF)  SAVED VECTOR OF LAGRANGE MULTIPLIERS.
*  RI  G(NF)  GRADIENT OF THE LAGRANGIAN FUNCTION.
*  RI  GO(NF)  SAVED GRADIENT OF THE LAGRANGIAN FUNCTION.
*  RO  R  VALUE OF THE STEPSIZE PARAMETER.
*  RO  F  NEW VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  OLD VALUE OF THE OBJECTIVE FUNCTION.
*  RO  P  NEW VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RI  PO  OLD VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RI  CMAX  VALUE OF THE CONSTRAINT VIOLATION.
*  RO  CMAXO  SAVED VALUE OF THE CONSTRAINT VIOLATION.
*  RO  DMAX  MAXIMUM RELATIVE DIFFERENCE OF VARIABLES.
*
* COMMON DATA :
*  II  NORMF  SCALING SPECIFICATION. NORMF=0-NO SCALING PERFORMED.
*         NORMF=1-SCALING FACTORS ARE DETERMINED AUTOMATICALLY.
*         NORMF=2-SCALING FACTORS ARE SUPPLIED BY USER.
*  II  ITERS  TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION.
*         ITERS=0 FOR ZERO STEP.
*
* SUBPROGRAMS USED :
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
*  S   MXVSET  INITIATION OF A VECTOR.
*  S   MXVSAV  DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE
*         SUBSTRACTED ONE.
*
      SUBROUTINE PYTRND(NF,N,X,XO,ICA,CG,CZ,G,GO,R,F,FO,P,PO,CMAX,CMAXO,
     & DMAX,KD,LD,ITERS)
      INTEGER NF,N,KD,LD,ITERS
      INTEGER ICA(*)
      DOUBLE PRECISION X(*),XO(*),CG(*),CZ(*),G(*),GO(*),R,F,FO,P,PO,
     & CMAX,CMAXO,DMAX
      INTEGER I,J,L
      DO 2 J=1,NF-N
      L=ICA(J)
      IF (L.GT.0) THEN
      CALL MXVDIR(NF,-CZ(J),CG((L-1)*NF+1),G,G)
      ELSE
      L=-L
      G(L)=G(L)-CZ(J)
      ENDIF
    2 CONTINUE
      IF (ITERS.GT.0) THEN
      CALL MXVDIF(NF,X,XO,XO)
      CALL MXVDIF(NF,G,GO,GO)
      PO=R*PO
      P=R*P
      ELSE
      F = FO
      P = PO
      CMAX = CMAXO
      CALL MXVSAV(NF,X,XO)
      CALL MXVSAV(NF,G,GO)
      LD = KD
      ENDIF
      DMAX = 0.0D0
      DO 1 I=1,NF
      DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D0))
    1 CONTINUE
      N=NF
      RETURN
      END