* SUBROUTINE PUDBG1                ALL SYSTEMS                92/12/01
* PORTABILITY : ALL SYSTEMS
* 92/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* VARIABLE METRIC UPDATE OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX
* USING THE FACTORIZATION B=L*D*TRANS(L).
*
* PARAMETERS :
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  RU  H(M)  FACTORIZATION B=L*D*TRANS(L) OF A POSITIVE
*         DEFINITE APPROXIMATION OF THE HESSIAN MATRIX.
*  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RA  S(NF)  AUXILIARY VECTOR.
*  RU  XO(NF)  VECTORS OF VARIABLES DIFFERENCE.
*  RI  GO(NF)  GRADIENTS DIFFERENCE.
*  RI  R  VALUE OF THE STEPSIZE PARAMETER.
*  RI  PO  OLD VALUE OF THE DIRECTIONAL DERIVATIVE.
*  II  NIT  ACTUAL NUMBER OF ITERATIONS.
*  II  KIT  NUMBER OF THE ITERATION AFTER LAST RESTART.
*  IO  ITERH  TERMINATION INDICATOR. ITERH<0-BAD DECOMPOSITION.
*         ITERH=0-SUCCESSFUL UPDATE. ITERH>0-NONPOSITIVE PARAMETERS.
*  II  MET1  SELECTION OF SELF SCALING. MET1=1-SELF SCALING SUPPRESSED.
*         MET1=2 INITIAL SELF SCALING.
*  II  MEC  CORRECTION IF THE NEGATIVE CURVATURE OCCURS.
*         MEC=1-CORRECTION SUPPRESSED. MEC=2-POWELL'S CORRECTION.
*
* SUBPROGRAMS USED :
*  S   MXDPGU  CORRECTION OF A DENSE SYMMETRIC POSITIVE DEFINITE
*         MATRIX IN THE FACTORED FORM B=L*D*TRANS(L).
*  S   MXDPGS  SCALING OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX
*         IN THE FACTORED FORM B=L*D*TRANS(L).
*  S   MXVDIF  DIFFERENCE OF TWO VECTORS.
*  RF  MXVDOT  DOT PRODUCT OF VECTORS.
*  S   MXVSCL  SCALING OF A VECTOR.
*
* METHOD :
* BFGS VARIABLE METRIC METHOD.
*
      SUBROUTINE PUDBG1(N,H,G,S,XO,GO,R,PO,NIT,KIT,ITERH,MET,MET1,MEC)
      DOUBLE PRECISION PO,R
      INTEGER          ITERH,KIT,MET,MET1,MEC,N,NIT
      DOUBLE PRECISION G(*),GO(*),H(*),S(*),XO(*)
      DOUBLE PRECISION A,B,C,GAM,PAR,DEN,DIS
      LOGICAL          L1,L3
      DOUBLE PRECISION MXVDOT,MXDPGP
      L1 = MET1 .GE. 3 .OR. MET1 .EQ. 2 .AND. NIT .EQ. KIT
      L3 = .NOT. L1
*
*     DETERMINATION OF THE PARAMETERS B, C
*
      B = MXVDOT(N,XO,GO)
      A = 0.0D0
      IF (L1) THEN
          CALL MXVCOP(N,GO,S)
          CALL MXDPGB(N,H,S,1)
          A=MXDPGP(N,H,S,S)
          IF (A.LE.0.0D0) THEN
              ITERH=1
              RETURN
          END IF
      END IF
      CALL MXVDIF(N,GO,G,S)
      CALL MXVSCL(N,R,S,S)
      C = -R*PO
      IF (C.LE.0.0D0) THEN
          ITERH = 3
          RETURN
      END IF
      IF (MEC.GT.1) THEN
          IF (B.LE.1.0D-4*C) THEN
*
*     POWELL'S CORRECTION
*
              DIS=(1.0D0-0.1D0)*C/(C-B)
              CALL MXVDIF(N,GO,S,GO)
              CALL MXVDIR(N,DIS,GO,S,GO)
              B=C+DIS*(B-C)
              IF (L1) A=C+2.0D0*(1.0D0-DIS)*(B-C)+DIS*DIS*(A-C)
          ENDIF
      ELSE
          IF (B.LE.1.0D-4*C) THEN
              ITERH = 2
              RETURN
          ENDIF
      ENDIF
      IF (L1) THEN
*
*     DETERMINATION OF THE PARAMETER GAM (SELF SCALING)
*
          IF (MET.EQ.1) THEN
              PAR = C/B
          ELSE IF (A.LE.0.0D 0) THEN
              PAR = C/B
          ELSE
              PAR=SQRT(C/A)
          ENDIF
          GAM = PAR
          IF (MET1.GT.1) THEN
              IF (NIT.NE.KIT) THEN
                  L3=GAM.LT.0.5D0.OR.GAM.GT.4.0D0
              ENDIF
          ENDIF
      ENDIF
      IF (L3) THEN
          GAM = 1.0D0
          PAR = GAM
      END IF
      IF (MET.EQ.1) THEN
*
*     BFGS UPDATE
*
      CALL MXDPGU(N,H,PAR/B,GO,XO)
      CALL MXDPGU(N,H,-1.0D0/C,S,XO)
      ELSE
*
*     HOSHINO UPDATE
*
      DEN=PAR*B+C
      DIS=0.5D0*B
      CALL MXVDIR(N,PAR,GO,S,S)
      CALL MXDPGU(N,H,PAR/DIS,GO,XO)
      CALL MXDPGU(N,H,-1.0D0/DEN,S,XO)
      ENDIF
      ITERH = 0
      IF (GAM.EQ.1.0D0) RETURN
      CALL MXDPGS(N,H,1.0D0/GAM)
      RETURN
      END